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Abstract. We consider the general problem of determining the steady state of 
stochastic nonequilibrium systems such as those that have been used to model (among 
other things) biological transport and traffic flow. We begin with a broad overview 
of this class of driven diffusive systems — which includes exclusion processes — focusing 
on interesting physical properties, such as shocks and phase transitions. We then 
turn our attention specifically to those models for which the exact distribution of 
microstates in the steady state can be expressed in a matrix product form. In addition 
to a gentle introduction to this matrix product approach, how it works and how it 
relates to similar constructions that arise in other physical contexts, we present a 
unified, pedagogical account of the various means by which the statistical mechanical 
calculations of macroscopic physical quantities are actually performed. We also review 
a number of more advanced topics, including nonequilibrium free energy functionals, 
the classification of exclusion processes involving multiple particle species, existence 
proofs of a matrix product state for a given model and more complicated variants of 
the matrix product state that allow various types of parallel dynamics to be handled. 
We conclude with a brief discussion of open problems for future research. 
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1. Introduction 



1.1. What is a nonequilibrium steady state? 

Consider a finite volume containing interacting particles and that is coupled at opposite 
boundaries to particle reservoirs held at different chemical potentials, as shown in 
Figure 1. In this situation one anticipates a net flux of particles from the reservoir 
with the greater chemical potential through the system and into the opposite reservoir. 
After some initial transients we expect a state to arise in which there is a nonzero 
mean flux that is constant over space and time. This flux reveals the system to be in 
a nonequilibrium steady state. More generally, we have in mind any physical system all 
of whose observables do not change with time, but nevertheless exhibits an irreversible 
exchange of heat, particles, volume or some other physical quantity with its environment. 

When heat, particle or volume exchanges are reversible, the system is at equilibrium 
with its environment and, as is well established, its microstates have a Gibbs-Boltzmann 
distribution. This knowledge allows one then to predict the macroscopic physics of a 
many-body system given a microscopic model couched in terms of energetic interactions. 
Furthermore, these Gibbs-Boltzmann statistics emerge from the application of a very 
simple principle, namely that of equal a priori probabilities for the system combined 
with its environment. 

By contrast, when a system is out of equilibrium, very little is known about the 
statistics of the microstates, not even in the steady state. This is despite a wide range of 
approaches aimed at addressing this deficiency in our understanding of nonequilibrium 
physics. These approaches fall roughly into two broad categories. 

First, one has macroscopic theories, such as Onsager and Machlup's pioneering 
work on near-equilibrium fluctuations [1,2]. In this regime, it is assumed that the 
thermodynamic forces that restore an equilibrium are linear, and it is possible to 
determine the probability of witnessing certain fluctuations in the equilibrium steady 
state, the most probable trajectory given by a principle of minimum energy dissipation. 
In this work, we are interested in steady states that are far from equilibrium, i.e., 
where the system is driven beyond the linear response regime. Recently there has been 
some success in extending the Onsager-Machlup theory of equilibrium fluctuations to 
such nonequilibrium steady states, as long as one has been able to derive macroscopic 
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Figure 1. A channel held out of equilibrium through the application of a chemical 
potential gradient. It is anticipated that a nonequilibrium steady state will be reached 
in which particles flow from one reservoir to another at a constant rate. 



5 



(hydrodynamic) equations of motion for the system [3]. Furthermore, over the past 
few years fluctuation theorems have been derived under a wide range of conditions. 
Typically, these relate the probabilities of entropy changes of equal magnitude but 
opposite sign occurring in a system driven arbitrarily far from equilibrium (see [4] for a 
brief introductory overview). Although some connections between these various topics 
have been established (see, e.g., [5]), it is fair to say that a complete coherent picture of 
the macroscopic properties of nonequilibrium steady states is still lacking. 

Somewhere between this macroscopic approach and a truly microscopic approach 
(i.e., one that would follow directly from some microscopic equations of motion) lie 
a range of mesoscopic models that are intended to capture the essential features of 
nonequilibrium dynamics — irreversibility, currents, dissipation of heat and so on — but 
are nevertheless simple enough that analytical treatment is possible. In particular, it 
is possible to explore the macroscopic consequences of the underlying dynamics, in the 
process identifying both similarities with equilibrium states of matter and the novel 
features peculiar to nonequilibrium systems. 

It is this approach that is the focus of this review article. Specifically, we will discuss 
a set of models that have a steady-state distribution of microstates that can be expressed 
mathematically in the form of a matrix product. We will explain in detail how these 
expressions can be used to calculate macroscopic steady-state properties exactly. These 
will include particle currents, density profiles, correlation functions, the distribution of 
macroscopic fluctuations and so on. These calculations will reveal phenomena that occur 
purely because of the far-from- equilibrium conditions: for example, boundary induced 
phase transitions and shock fronts. Meanwhile, we will also find conceptual connections 
with equilibrium statistical physics: in some of the models we discuss, for example, 
it is possible to construct partition functions and free-energy-like quantities that are 
meaningful both for equilibrium and nonequilibrium systems. We begin by explaining 
this modelling approach in more detail. 

1.2. Modelling nonequilibrium dynamics 

In Figure 2 we have sketched some molecular motors (kinesins) that are attached to a 
microtubule. This is essentially a track along which the motors can walk by using packets 
of energy carried by ATP molecules present in the surroundings. Clearly, at the most 
microscopic level, there are a number chemical processes at play that combine to give 
rise to the motor's progress along the track. At a more coarse-grained, mesoscopic level, 
we can simply observe that the motor takes steps of a well-defined size (approximately 
8nm for a kincsin [6]) at a time, and thus model the system as a set of particles that hop 
stochastically between sites of a lattice. This stochastic prescription is in part intended 
to reflect the fact that the internal degrees of freedom that govern when a particle hops 
are not explicitly included in the model. 

The simplest such stochastic model for the particle hops is a Poisson process that 
occurs at some prescribed rate. Let C and C be two conflgurations of the lattice that 
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Figure 2. Kinesins attached to a microtubule that move in a preferred direction by 
extracting chemical energy from the environment. At a mesoscopic level, this can be 
modelled as a stochastic process in which particles hop along a one-dimensional lattice 
as shown. 

differ by a single particle hop. We then define W{C — > C) as the rate at which this hop 
occurs, such that in an infinitesimal time interval dt the probability that that hop takes 
place is W{C — >• C')dt. The rate of change of the probability P{C,t) for the system to 
be in configuration C at time t is then the solution of the master equation 

^P(C, t) = ^ P(C', t)W{C' - C) - 5^ P(C, t)W{C ^ C) (1.1) 
CYC cvc 

subject to some initial condition P(C,0). The first term on the right-hand side of this 
equation gives contributions from all possible hops (transitions) into the configuration 
C from other configurations C; the second term gives contributions from hops out of C 
into other configuration C. We are concerned here with steady states, where these gain 
and loss terms exactly balance causing the time derivative to vanish. Such solutions of 
the master equation we denote P(C). In this work, nearly all the models we discuss have 
the property that a single steady state is reached from any initial condition: i.e., they 
are ergodic [7,8]. Given then that the distribution P(C) is unique for these models, we 
concentrate almost exclusively on how it is obtained analytically, and how expectation 
values of observables are calculated. The remaining models that are not ergodic over 
the full space of configurations turn out to have a unique steady-state distribution over 
some subspace of configurations. We will see that such a distributions of this type can 
also be handled using the methods described in this work. 

1.3. A catalogue of steady states 

Solutions of the master equation (1.1) can assume a number of structures. Some general 
classes of steady state — not intended to be mutually exclusive — can be summarised as 
follows: 



Equilibrium (Gibbs-Boltzmann) steady state If the internal energy of a microstate C 
is E{C) and the system as a whole is in equilibrium with a heat bath with inverse 
temperature /3, the microstates have a Gibbs-Boltzmann distribution P(C) oc e"'^^''''^ 
The dynamics is reversible in that the probability of witnessing any particular trajectory 
through phase space is equal to that of its time reversal. This then implies that stochastic 
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dynamics expressed in terms of transition rates — > C) must satisfy the detailed 
balance condition [9, 10] 

P{C)W{C C) = P{C')W{C' C) (1.2) 

where here P{C) is the Gibbs-Boltzmann distribution. A simple consequence of (1.2) is 
that there are no fluxes in the steady state. Hence any system that does exhibit currents 
in the steady state, i.e., nonequilibrium systems, generally do not have a stationary 
distribution that satisfies detailed balance. 

Factorised steady state Sometimes, whether in or out of equilibrium, one has a 
factorised steady state. Typically one has in mind a lattice (or graph) with sites with 
configurations specified by occupancy variables rii (i.e., site i contains n,j particles). A 
factorised steady state then takes the form P{C) oc YliLi fii'i^i)- such a structure 
to arise, certain constraints on the transition rates W{C C) must be satisfied. 
For example, in an equilibrium system, one would need the energy of the system to 
be the sum of single-site energies. Out of equilibrium, one finds factorised steady 
states in zero-range processes, in which particles hop at a rate that depends only 
on the occupation number at the departure site. Moreover, necessary and sufficient 
conditions for factorisation in a broad class of models which includes the ZRP have been 
established [11]. Zero-range processes and various generalisations have been reviewed 
recently in a companion paper [12] and so we do not discuss them further here. 

Matrix product steady state A matrix product steady state is an extension to a 
factorised steady state that is of particular utility for one-dimensional models. The 
rough idea is to replace the scalar factors fi{ni) with matrices the steady state 
probability then being given by an element of the resulting matrix product Ilili^ni- 
Since the matrices Xn for different occupancy numbers n need not commute, one opens 
up the possibility for correlations between the occupancy of different sites (above those 
that emerge from global constraints, such as a fixed particle number, that one sees in 
factorised states). It turns out that quite a number of nonequilibrium models have a 
matrix product steady state, and we will encounter them all in the course of this review. 

Most general steady state In principle, the (assumed unique) stationary solution of the 
master equation (1.1) can always be found if the number of configurations is finite, since 
(1.1) is a system of linear equations in the probabilities P{C). One way to express this 
solution is in terms of statistical weights f{C), each of which is given by the determinant 
of the matrix of transition rates obtained by removing the row and column corresponding 
to the configuration C (see e.g. [13,14]). To arrive at the probability distribution P{C) 
one requires the normalisation Z = X]c/(^)' ^'^ th.at then P{C) = f{C)/Z. This 
normalisation has some interesting properties: it is uniquely defined for any ergodic 
Markov process with a finite number of configurations; it can be shown to be equal to the 
product of all nonzero eigenvalues of the matrix of transition rates; and it is a partition 
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function of a set of trees on the graph of transitions between different microscopic 
configurations. In this latter case, the densities of particular edges in the ensemble of 
trees are controlled by tuning the transition rates corresponding to those edges, so in this 
interpretation these transition rates are equivalent to equilibrium fugacities. Although 
this interpretation of the normalisation is rather abstract, concrete connections between 
transition rates in certain nonequilibrium models that have a matrix product steady 
state and fugacities in an equilibrium ensemble have been made, and these we shall 
discuss later in this review. We also remark that the equilibrium theory of partition 
function zeros to characterise phase transitions also carries over to the normalisation as 
just defined, at least for those models that have been tested [15-17]. 

1.4- Purpose of this review 

We have three principal aims in this work. First, we wish to illustrate the insights 
into nonequilibrium statistical mechanics that have been gained from exactly solvable 
models, particularly those that have a steady state of the matrix product form. Secondly, 
we seek to provide a self-contained pedagogical account of the various analytical methods 
and calculational tools that can be used to go from the matrix product expressions to 
predictions for the macroscopic physics. Finally, we wish to present a thorough review 
of the progress that has been made using the matrix product approach over the last 
few years. As far as we are aware, the matrix product approach has not been the 
focus of a review for nearly a decade [18, 19], and we feel that it is high time that 
the significant developments that have occurred in the meantime should be collected 
together in one place. In order to prevent this review from becoming unmanageably 
long, we focus purely on static physical properties exhibited in the steady states of 
models solvable by the matrix product method. This means we have unfortunately had 
to omit discussion of some very interesting topics, for example how some dynamical 
properties, such as fluctuation phenomena, have been elucidated through the use of the 
Bethe ansatz [20-22], determinental solutions and their connection to random matrix 
theory [23-25] and dynamical matrix products [26,27]. Some of these topics, however, 
have recently been reviewed elsewhere [28,29]. We also direct the reader to established 
reviews, such as [19,30-35], for any other background that we have been forced to leave 
out here. 

To realise the aims stated above, we provide in the next section a general account 
of the physics one expects to see in nonequilibrium dynamical systems, and outline the 
essential ideas underlying the matrix product approach. Thereafter, we go into the 
details of how the simplest models are solved, and show a number of contrasting (but 
equivalent) approaches that have found apphcation to more complex problems. These 
cases then form the material of the remainder of the review, which we round off by 
posing some open problems for future research. 
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2. An overview of driven diffusive systems 

2.1. A paradigmatic example: the asymmetric simple exclusion process (ASEP) 

The asymmetric simple exclusion process (ASEP) is a very simple model of a driven 
system in one dimension that has biased diffusion of hard-core particles. We shall discuss 
two versions: the periodic system and the open system. The latter is coupled to particle 
reservoirs at either end so that, as described in the previous section, there is a steady 
state with a constant particle flux. The open system has a long history, having first 
appeared in the literature to our knowledge as a model of biopolymerisation [36] and 
transport across membranes [37]. In the mathematical literature meanwhile diffusion 
with collisions between particles was first studied by Harris [38] and the terminology 
simple exclusion was first defined by Spitzer [39]. Over the years, applications to other 
transport processes have appeared, e.g., as a general model for traffic flow [40] and 
various other theoretical and experimental studies of biophysical transport [41-45]. 
Whilst interesting and important, these many applications arc not our primary concern 
here. Also we shall not do justice to the many rigorous mathematical results which 
have been summarised in the books by Liggett [46,47]. Rather, our interest in the 
ASEP lies in its having acquired the status of a fundamental model of nonequilibrium 
statistical physics in its own right in much the same way that the Ising model has become 
a paradigm for equilibrium critical phenomena. In particular, the ASEP and related 
models show — despite their supcrflcial simplicity — a range of nontrivial macroscopic 
phenomena, such as phase transitions, spontaneous symmetry breaking, shock fronts, 
condensation and jamming. 

In common with the vast majority of the models we discuss in this review, the 
ASEP is deflned as a stochastic process taking place in continuous time on a discrete 
one-dimensional lattice of N sites. Hopping is totally asymmetric: a particle sitting 
to the left of an occupied site hops forwards one site (as in Figure 2), each hop being 
a Poisson process with rate p. That is, in an inflnitesimal time interval di, there is a 
probability pdt that one of the particles that can hop, does hop, the identity of that 
particle being chosen at random. In simulation terms, this type of dynamics corresponds 
to a random- sequential updating scheme, in which bonds between particles are chosen 
at random and then, if there is a particle at the left-hand end of the bond and a vacancy 
at the right, the particle is moved forwards. In this scheme, each update corresponds 
(on average) to l/{Np) units of time; in the following we take p — 1 without loss 
of generahty. The deflnition of the model is completed by specifying the boundary 
conditions. 

2.1.1. Periodic boundary conditions The simplest version of the model has a ring of N 
sites as depicted in Figure 3(a): a particle hopping to the right on site lands on site 
1 if the receiving site is empty. Under these dynamics, the total particle number M is 
conserved. Furthermore, as we will show in Section 3.1, the steady state is very simple: 
all conflgurations (with the allowed number of particles) are equally likely. As there 
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Figure 3. Asymmetric exclusion process with (a) periodic and (b) open boundary 
conditions. The labels indicate the rates at which the various particle moves can 
occur. 
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are (^) allowed configurations the probability of any one is The steady-state 

current of particles, J, through a bond i + 1 is given by p (here equal to 1) multiplied 
by the probability that there is a particle at site i and site i + 1 is vacant. One finds 

(mIi) _ M (N - M ) 

in the thermodynamic limit where A^, M — oo at a constant ratio M/N = p. 

One may also ask, starting from a known configuration, how long the system takes 
to relax to the steady state. It turns out that this timescale T depends on the system 
size via the relation T ~ A^^ which defines a dynamic exponent z. One can show 
using the Bethe ansatz (not covered here, but see, e.g., [20,48,49]) that for the ASEP 
on a ring, z = 3/2. Also in the symmetric case where particles hop both to the left and 
to the right with equal rates, but still with exclusion, one obtains z = 2 which is the 
value one would find for a purely diffusive process. 



2.1.2. Open boundary conditions To model the interaction of an open system with 
reservoirs at different densities, Poisson processes acting at the boundary sites are added. 
Specifically, a particle may be inserted onto the leftmost site with rate a (if it is vacant), 
and may leave the system from the rightmost site at rate (3: see Figure 3(b). Here we 
will consider only a, (3 < 1 although there is no problem in considering rates outwith 
this range. Then we may think of reservoirs at site with density a and a reservoir at 
site A^+1 with density l — (3. To include these moves in the simulation scheme described 
above, one would additionally allow the bonds between the system and the reservoirs 
to be chosen, and perform the particle updates with probability a (entry) or (3 (exit) 
according to the bond chosen (assuming a and (3 both less than 1). These moves admit 
the possibility of a nonzero current in the steady state. 

The open system is more interesting in that phase transitions may occur in the 
steady state. In the following we shall explore the origins of these phase transitions, both 
through approximate treatments and an exact solution. In this context we recognise — 
by direct analogy with equilibrium phase transitions — a phase transition as a sudden 
change in form of macroscopic quantities such as the particle current across a bond. 
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or the density at a site. It turns out that the particle current plays a central role, 
analogous to the free energy of an equilibrium system, in determining the nature of the 
phase transition i.e. the order of the transition is determined by which derivative of the 
current with respect to some external parameter exhibits a discontinuity. 

We remark that other types of updating scheme can also be considered. For 
example, in applications to traffic flow or pedestrian dynamics [33,50], it is more 
natural to consider parallel dynamics where many particles can hop in concert. The 
ASEP remains solvable under certain classes of parallel updating schemes which will be 
discussed towards the end of this review (Section 10). 



2.2. Mean-field and hydrodynamic approaches 

2.2.1. Lighthill-Whitham theory of kinematic waves We begin by reviewing a classical 
phenomenological theory, first appfied to traffic flow [51], which serves as a flrst recourse 
in our understanding of the phase diagrams of driven diffusive systems. The idea is to 
postulate a continuity equation for the local density p 

dp dJ ^ . . 

where J is the current of particles. (Note we use here a different nomenclature and 
notation to [51] which discusses concentration and ffow rather density and current.) 

The analysis rests on the key assumption that there is a unique relation, J(p), 
between the current and local density. Also, it is assumed that there is a maximum 
current Jm (the capacity of the road) at density pm- The flrst assumption implies that 
(2.2) becomes 

(-) 

where 

<P) - ^ . (2.4) 

Now, an implicit solution of (2.3) is of travelling- wave form 

p{x,t)^f{x-Vg{p)t) (2.5) 

as can readily be checked by substitution into (2.3). The arbitrary function / is 
determined by the initial density proflle p(x, 0) = fi^)- The interpretation of the 
solution is that a patch with local density p propagates with velocity Vg{p). The 
propagation of such a patch is known as a collective phenomenon known as a kinematic 
wave and Vg is a group velocity. The velocity c of a single particle in an environment of 
density p, on the other hand, is deflned through 

J(p) = c{p)p . (2.6) 

One sees that 

d , , dc ,^ ^, 

., = -(pc) = c + p-. (2.7) 
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Figure 4. Kinematic waves and sharpening of a shock: the figure sketches the 
evolution of a density profile p{t) from an initial profile p(0) according to (2.2) and 
illustrates how sharp discontinuities in the profile may develop (see text for discussion). 



Thus, if the single particle velocity is a decreasing function of density, which is what we 
expect, we find Vg < c and kinematic waves travel backwards in the frame of a moving 
particle. 

If Vg is a decreasing function of p we have the phenomenon of shock formation (see 
Figure 4). That is, since the patches of an initial density profile travel with different 
speeds, the low-density regions catch up with higher density regions and discontinuities 
in the density profile, known as shocks, develop. Strictly, after the formation of a shock 
the description of the density profile by the first order equation (2.2) breaks downs and 
one has to supplement (2.2) with second-order spatial derivatives to describe the shock 
profile. 

To deduce the velocity of a shock consider two regions of density pi to the left 
and p2 to the right separated by a shock. Mass conservation implies that the velocity 
(positive to the right) of the shock is given by 

V, = ^(^^^ - ^^''^ . (2.8) 

Note that if the current density relation has a maximum then one can have a stationary 
shock, Vs = 0. 

A particular form for the current density relationship studied in [51] and which is 
relevant to the ASEP is 

J(p)=p(l-p). (2.9) 

In this case = 1/4 at = 1/2 and J is symmetric about the maximum. Expression 
(2.9) may be easily understood for the ASEP by noting that in order for a particle to 
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hop across a bond, the site to the left must be occupied and that to the right empty. 
For a bond at position x, the probabihty of the former is p{x) and that of the latter 
is 1 — p{x). Assuming these events are uncorrelated (which they are not — but we shall 
come to this later), we have J(p) = p(l — p). 

In the case (2.9) one can compute the kinematic wave and shock velocities from 
(2.4,2.8) 

Vg = l-2p Vs = l-Pi-P2 (2.10) 

where, as above, p\ and p2 are the densities on either side of the shock. Thus the 
kinematic wave velocity is negative when p > 1/2 and the shock velocity is negative 
when pi > 1 — p2. 

The kinematic wave theory can be used to predict the phase diagram of the open 
boundary ASEP, in which distinct steady-state behaviours are demarcated (Figure 5). 
The left hand boundary x = is considered as a reservoir of particle of density pi — a 
and the right hand boundary {x — N) as a, reservoir of density Pr — 1 — (3- Associated 
with these boundary densities are kinematic waves with velocities 

vi^l-2a ; Vr^2p-1 . (2.11) 

In the case where a < 1/2 and P < 1/2 both kinematic waves propagate into the system. 
So, for example, from an initially empty system the kinematic waves of densities a and 
1 — /3 will enter from the left and right of the system and meet somewhere in the middle 
forming a shock which then moves with velocity 

v, = (3-a. (2.12) 

li (3 > a the shock moves to the right hand boundary and density associated with the 
left hand boundary, p; = a, is adopted throughout the bulk of the system. On the other 
hand if a < /3 the shock moves to the left hand boundary and density associated with 
the right hand boundary, pr — 1 — P, is adopted throughout the bulk of the system. 

In the case a — P < l/2,Vs — and the shock is stationary. In the stochastic system 
the shock, although on average stationary, diffuses around the system and effectively 
reflects off the boundaries. The result is that the shock is equally hkely to be anywhere 
in the system. 

In the case where one of a or > 1/2, the kinematic wave associated with that 
boundary does not propagate into the system and the kinematic wave which does 
propagate from the other boundary controls the bulk density. Thus, the boundary 
with a or P < 1/2 controls the bulk density 

Finally if both a, P > 1/2. The kinematic waves from both boundaries do not 
penetrate. To describe this phase one needs to add a diffusive contribution to the 
current (2.9) i.e. to consider second-order spatial derivative of p which we shall do in 
the next section. The result is that steady state of the system has density 1/2; the 
system adopts the maximal current density pm — 1/2 which is the density associated 
with kinematic wave velocity zero. 
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Region 


Phase 


Current J 


Bulk density p 


a < (3,a < ^ 


Low-density (LD) 


a{l — a) 


a 


/3 < < i 


High- density (HD) 


(3{l-(3) 


l-(3 


a > 1,(3 > I 


Maximal current (MC) 


1 

4 


1 

2 



Table 1. Properties of the ASEP obtained using the extremal-current principle under 
a mean-field approximation for the density-current relationship. 




Figure 5. Mean field phase diagram for the ASEP showing the range of a and (3 in 
which the low-density (LD), high-density (HD) and maximal-current (MC) phases are 
seen. 



The resulting bulk densities and currents are then as shown in Table 1 which 
corresponds to the phase diagram given in Figure 5. Here we have adopted what 
has become the standard nomenclature for the three phases — the terminology should 
hopefully be self-explanatory. 

This simple kinematic wave theory using the first order equation (2.2) correctly 
predicts the phase diagram, as we shall see from the exact solution. The theory also 
implies two sub-phases of the high and low-density phases according to whether the 
boundary which does not control the bulk density has a or /3 greater than or equal to 
1/2. 

A more refined theory is the domain wall theory [52-55] which is able to make 
more general predictions (such as for multi-species problems) and is reviewed in [35]. 
We remark that kinematic waves are in fact the characteristic curves of the differential 
equation (2.3). For more general first-order differential equations modelling driven 
diffusive systems one can use the method of characteristics to predict the phase diagram 
as we briefly discuss in Appendix A. 

2.2.2. Comparison of mean-field and hydrodynamic approaches In the previous 
subsection we showed how a simple phenomenological description by a first-order 
differential equation can predict the correct phase diagram of the ASEP. Generally the 
assumption that the current is solely a function of the local density is an assumption 
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which is mean-field in nature. Mean-field theories are usually considered as a general 
class of approximations where correlations within the configuration of the system 
are ignored at some level. Thus a very simple mean-field theory is to only keep 
information about the density profile. A more sophisticated mean-field theory might 
keep information about two point correlation functions but will still ignore correlations 
at some level. In some cases or limits a mean-field description may become an exact 
description. 

At this point, it is useful to compare the idea of a mean- field description with the 
hydrodynamic limit usually considered in the mathematical literature. In the latter 
case one defines, for example, a local density by averaging or coarse graining over some 
scale. Then one seeks to define a limit where this scale is allowed to become large 
but the ratio of this scale to the system size is small. The hydrodynamic limit is then 
proven by showing that the variables coarse-grained on this scale satisfy a deterministic 
equation where the noise has been scaled out. The procedure involves determining the 
local stationarity of the system on the coarse-grained scale. As in mean-field theory, 
the simplest case results in describing the system by the density profile alone without 
correlations. 

Therefore, at a practical level, the results of a mean-field description often coincide 
with those obtained in a hydrodynamic limit — both approaches result in the same 
differential equations involving the density or some small number of variables. But, it is 
important to stress that the approaches differ significantly in spirit: the hydrodynamic 
limit specifies a limit or scale where an exact equation is proven. On the other hand, 
mean-field theory is often used as a rough and ready approximation regardless of whether 
the hydrodynamic limit exists. In addition, the hydrodynamic limit introduces coarse 
grained densities whereas the equivalent mean-field equations are often obtained without 
such an explicit procedure as we show below. Finally mean-field or hydrodynamic 
equations can often be written down phenomenologically as with the Lighthill-Witham 
theory. 

2.2.3. Mean-field treatment of the density profile We start our considerations of mean- 
field approximations by keeping space discrete, and implementing a simple procedure 
which is to to assume a factorised form for the stationary distribution. That is, in terms 
of the indicator variables Tj which take values Tj = 1 if site i is occupied and Ti — 
otherwise, we suppose 



and Pi = (Ti). Note that all higher order correlations vanish i.e. {riTj) is approximated 
in this mean-field theory by piPj. 

Using the master equation (1.1), one can derive an exact equation for the evolution 
of the density at site i: 




i=l 



N 



where 





(2.14) 
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Figure 6. Iteration of the map (2.18) in the high and low-density phases. Profiles 

A and B iterate towards the high-density fixed point; profiles C and D begin 
infinitesimally close to the low-density fixed point and iterate away from it as the 
right boundary is approached. 



A full derivation is given, for example, in [18]. To understand this equation note that 
the right-hand side contains the particle current from site i — 1 to i minus the particle 
current from itoi + 1, which by particle conservation gives the rate of change of density 
at site i. In the steady state, where the density is independent of time, the right hand 
side of (2.14) must vanish which implies the exact steady state result 

J = {nil - n+i)) (2.15) 

where the current J between two neighbouring sites is independent of position. 
In the mean-field approximation equation (2.14) becomes 



pi) - pi{l - pi+i) (2.16) 
and the steady-state mean-field current becomes 

J = p,(l - p,+i) . (2.17) 

The density profile and phase diagram can be obtained in the mean-field approximation 
by considering the mapping 

J 



1 - 



The mapping has fixed points at 
1 
2 



P± 



1±\/1-4J 



(2.18) 



(2.19) 



The structure of possible solutions is illustrated in Figures 6 and 7 and correspond to 

the three phases. In the high-density phase the map iterates to a value of pi arbitrarily 
close to the stable p_|_ fixed point in some finite number of steps. In the low-density 
phase the map starts close to the unstable p_ fixed point. In the maximal current phase 
J = 1/4 + 0{N~'^) so that there arc no real fixed points and the density profile takes 
on its characteristic maximal current shape: the gradient of the density profile is always 
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Figure 7. Iteration of the map (2.18) to give the density profile in the maximal-current 
phase. 



negative and far away from both boundaries the density approaches the value 1/2 (see 
Figure 7). 

The phase diagram can be fully analysed within this discrete mean-field 
approximation [56]. Here we choose to make a further continuum approximation as 
this approach can be easily generalised to other systems as we shall see in the next 
subsection. One can take a continuum limit of (2.16) in a simple way by replacing pi{t) 
by p{x, t) where x — ai and a is introduced as the lattice spacing. Expanding to second 
order in a 



yields 



where 



Pi±i 
dp 

at ^ 



_^ dp d^p 



—a 



d_ 

dx 



Up) - D 



dp 
dx 



Jo{p) = p{l - p) and D = 



(2.20) 



(2.21) 



(2.22) 



Jq{p) represents a contribution to the current due to the asymmetry of the dynamics i.e. 
the external drive and —D^p{x) is the diffusive contribution to the current (down the 
concentration gradient). An equation of this form was postulated by Krug [57]. Ignoring 
the diffusive part recovers the Lighthill-Whitham phcnomcnological equation (2.2) or 
more formally rescaling time to t' = t/a and letting a ^ recovers the hydrodynamic 
limit on the Euler scale. Retaining the diffusive term in (2.21) allows one to calculate 
the (approximate) density profiles explicitly. To see this one considers the steady state 
of (2.21) 



J = J^(p) - (2.23) 
where J is a constant (the current) to be determined. Rearranging (2.23,2.22) yields 

(2.24) 
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dp 
dx 



{p - P+){p - P-) 
D 




Figure 8. Sketches of the possible solutions of equation (2.24) with p± real, which 
corresponds to J less than the maximum of Jo{p). 



where p± are given by (2.19). One can separate variables of this first order differential 
equation and integrate, for example, from the left boundary 



p{x) 



dp 



D 



{P - P+)(P - P-) 



— X 



which yields 



In 



{p-p+){piO)-p-) 
(p-p_)(p(0)-p+) 



D 



(2.25) 



(2.26) 



The boundary conditions to be imposed on (2.23) are implied by the reservoir densities: 
p(0)=a p{N + 1) = 1-P (2.27) 
the second of these will fix J from (2.26): 



{l-f3-p+){a-p.) 
(l-/5-p_)(a-p+)J 



= exp 



(p+-p-)(iV+l) 
D 



(2.28) 



One can rearrange (2.26) to give an explicit expression for the profile, but to deduce 
the phase diagram it is simplest to return to (2.24). One can easily sketch the possible 
solution curves to this equation as is done in Figures 8 and 9. The phase diagram in 
the hmit of a large system can be deduced as follows: 

High-density phase: The bulk density is p+ = p{N-\-l) = 1 — /? which yields J = /3(1 — /3) 
and p_ = /9. Then the condition p+ > p_ implies (3 < 1/2 and the condition that 
p(0) > p_ yields a < (3. 

Low-density phase: The bulk density is p_ = p(0) = a which yields J = a{l — a) and 
p+ = q;(1 — a). Then the condition p+ > p_ implies a < 1/2 and the condition that 
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Figure 9. Sketches of the solutions of equation (2.24) with p± complex, which 
corresponds to J greater than the maximum of Jo{p)- 



p{N + 1) < p+ yields a > p. 



Co-existence line: When a = f3 < 1/2 the profile begins infinitesimally close to p- — a 
and finishes infinitesimally close to p+ = 1 — /9. In the bulk of the system there is a 
shock front in which the density switches over. Thus, on this line the High-density and 
Low-density phases coexist. 

Maximal current phase: In this case J > 1/4 and the profile has the shape shown in 
Figure 9. One can show that to satisfy (2.28) one must take 

Then setting p{x) — l/2 + r]{x), (2.24) becomes 

| = 4-0(l/7V^) (2.30) 

which implies that in the large N ^ oo limit the decay of the density away from the 
left boundary follows, for large x, 

p{x) - 2 ~ - ■ (2.31) 

2.2.4- Extremal current principle Let us consider more general equations of the form 
(2.23) where the current of particles J across each bond is the sTim of a systematic 
current Jo{x) and a diffusive current Jo- The former is taken to be a known function 
Jo(p) of the local density p{x) arising from the nature of the particle interactions whilst 
the latter is taken to cause motion down a density gradient. Thus the diffusive current 
should depend only on the local gradient, take opposite sign to its argument and vanish 
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in regions of constant density: e.g. Jo = —D{p)p' where D{p) is some (positive) diffusion 
constant. 

The steady-state profile is given by (2.24) where J is a constant and D may now 
depend on p. Generally, if Jq has a single maximum J™"^ as a function of p we obtain 
density profiles with the same qualitative structure as just deduced for the ASEP. To 
see this we can again sketch the solution curves as in Figs. 8 and 9. If J < J^"^ there 
are two fixed points p± where = and the high and low-density profiles iterate 
to and away from these fixed points. If J > J^"^ there are no real fixed points and 

= yielding a maximal current profile. As the boundary reservoir densities a and 
1 — /9 are varied, the system undergoes transitions, of the same type as those exhibited 
by the ASEP, between regimes characterised by different functional forms of the particle 
current and density. 

These phase transitions can be understood qualitatively using a macroscopic 
extremal current principle proposed by Krug [57] and later developed into a more general 
description of driven systems [53,54]. The argument runs as follows. Prom Figs. 8 and 
9 one sees that the density profile has the following properties: it is monotonic in x] 
it is bounded from above and below due to the finite values of the boundary reservoir 
densities, and in the limit of an infinite system there is a bulk region in which the density 
is roughly constant p{x) = p. In the bulk region, the diffusive current vanishes, and 
so here J = Jo{p)- If the boundary conditions are such that the density at the left 
boundary is greater than that at the right Pr, we have regions near the boundaries 
where the density is a decreasing function of x (as a consequence of monotonicity) . In 
these regions, the diffusive current Jd{p) is positive, and hence the systematic current 
is reduced compared to that in the bulk. This implies that 

J = Jo{p) = max Jq{p) when pn < pL (2.32) 

since every possible density p in the range [pR, pi\ is represented somewhere in the 
system, and as we just argued, the systematic current in the bulk must be greater than 
that at any other point. A similar argument goes through for the case pl < Pr, except 
here the diffusive current is negative in the boundary regions and hence 

J = Jo(p) = min Jo(p) when pi < Pr ■ (2.33) 

The pair of equations (2.32) and (2.33) encapsulate the extremal current principle. 

2.3. The exact solution of open-boundary ASEP 

So far we have discussed mean-field approaches to deducing the phase diagram of 
the ASEP. These approaches yield the correct phase diagram but do not correctly 
predict correlation functions. For example, the density profiles discussed in the previous 
subsections are only qualitatively correct. The time has now come to summarise 
the main ideas behind the exact solution of the ASEP that uses the matrix product 
formulation briefiy introduced in Section 1 [58]. (We note that an exact solution using 
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Figure 10. The procedure by which a particle configuration is transhterated to a 
string of matrices. 



recursion relations for the stationary weights was found independently [59] , building on 
earlier work [56,60] for the case a = f3 = 1.) 

2.3.1. Nonequilihrium matrix product steady states To recap, a matrix product state 
is constructed as a product of matrices, one for each site, chosen according to the 
state of the site. Then, a scalar probability is obtained by combining in some way the 
elements of this product. The precise prescription for the ASEP with open boundaries 
was determined in [58]. It is summarised as follows. 

First, it is convenient to work with an unnormalised weights /(ri, . . . , tat) so that 

P(ri,r2,...,r7v) = (2.34) 

in which is the normalisation obtained by summing the weights / over all 2^ possible 
configurations of the A^-site lattice. It is these weights that take the matrix-product form 

/^(n, r2, . . . , rjv) = mXr,X^, ■■■XrAV) (2.35) 

in which X^-. is a matrix D if site i is occupied by a particle (x; = 1) or a matrix 
E otherwise. Observe that this is a very visual way of representing a particle 
configuration — see Figure 10. Note further that the matrix appearing at a particular 
point in the product does not depend explicitly on the site label only the state of that 
site, and that the vectors {W\ and \V) perform the necessary conversion of the matrix 
product to a scalar. 

With the notation established, expressions for macroscopic quantities quickly follow. 
In all of these, one requires the normalisation which can be written as 

Zn = {W\{D + E)''\V) = {W\C''\V) (2.36) 

in which the matrix C = D + E has been introduced for convenience. The quantity Zn 
is simply the sum of the weights of all possible configurations. The density at site i is 
then 

{W\C^-'DC^-^\V) 
Pi = = ^ (^•'^'J 

Zjn 

and the current between bonds i and i + 1 is 

Ji,i+i = {ri{l - Ti+i)} = . (2.38) 

Zjn 
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In order actually to evaluate these expressions, one needs to know a set of algebraic 
relations between D, E, \V) and {W\. In [58] it was shown that sufficient conditions for 
the stationary solution of the master equation of the ASEP via (2.35) are 

DE ^D + E (2.39) 

D\V) = ^\V) (2.40) 

{W\E^-{W\. (2.41) 
a 

In Section 3.1 we will prove that these relations indeed give the desired solution of 
the master equation. In the meantime, we can at least check the necessary condition 
that the steady state particle current Ji,i+i across the i*^ bond is independent of i, and 
further equals the currents of particles into and out of the system at its boundaries. To 
see this note first of all that the combination DE appearing in (2.38) can be replaced 
with D + E = C using (2.39). Then, 

At the left boundary we have 

Jl = ail - pi) = ' ^ ^ = ^ ' ' ^ = J (2.43) 

Zn Zm 

where we have used (2.41). Similarly, at the right boundary one has, using (2.40) 

{W\C^-'D\V) {W\C^-'\V) 

Jr = PPn = P ^ = ^ = J ■ (2.44) 

Zn Zn 

Relations (2.39)-(2.41) form what is referred to as a quadratic algebra which 
amounts to a set of reduction relations for matrix products. The power of the quadratic 
algebra is that it allows one to consider any arbitrary product U oi D and E matrices 
and determine the element (VF|C/|y) without recourse to an explicit representation. 
This is because if one has a product in which a D matrix appears before an E matrix, 
it will be possible to use (2.39) to reduce to the sum of two shorter products, one of 
which is missing the D and the other the E. If either of these terms themselves has a 
D before an E, (2.39) can again be used to generate further terms in the sum. Once 
all terms contain only D matrices appearing after E matrices, or completely lack either 
Ds or Es no further reduction is possible and one has found an expression for U in the 
form 

C/ = J^an.m^'" (2.45) 

n,m 

where a„ m is a positive coefficient arising from the reduction process. Once this has 
been achieved, evaluating the element using (2.40) and (2.41) is simphcity 

itself: 



{W\U\V) = Va„,„-— {W\V) . (2.46) 
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Conventionally, the combination is taken as unity, although its actual value never 

enters into any physical quantities, since it appears as a prefactor in both the numerator 
and denominator of (2.34). 

One might ask whether the order with which (2.39) is applied to different DE pairs 
appearing in U could affect the final outcome (i.e., whether (2.39)-(2.41) describes a 
non-associative algebra). That this is not the case can be seen from a straightforward 
argument. Let U comprise at least two DE pairs: U — Ui{DE)U2{DE)Us where C/1,2,3 
are products of D and E matrices that may, or may not, contain DE pairs. Clearly, 
the same matrix product expression is obtained by applying (2.39) to the two DE pairs 
in either order. Therefore, any non-uniqucncss would have to come at the next stage 
where there are four strings to reduce. Now, each of these strings contains either no 
DE pairs (which is unique); one DE pair (which can be reduced only in one way); or 
is in the form given for U above. In the latter case, one can iterate the argument just 
given to demonstrate that the reduction is unique. Hence, the quadratic algebra (2.39)- 
(2.41) involving only two matrices D and E is associative. However, as we shall see in 
Section 7, with more than two matrices the requirement for uniqueness under reduction 
is important and imposes restrictions on the class of models that can be solved using 
the matrix product method. 

The unconvinced reader may take refuge in the fact that explicit representations of 
the D and E matrices can be constructed — for example, that presented in Section 3.2.2. 
In general, these matrices do not commute and have infinite dimension, as was shown 
in [58]. An exception to this can be determined by assuming that D and E do commute: 
then, using (2.39)-(2.41) one finds that 



which can be true only if a + /5 = 1. Along this line, the relations (2.39)-(2.41) can be 
satisfied by the scalar choices D — 1/P,E = l/a, {W\ = \V) — 1. This is equivalent to 
putting Pi — a/{a + 13) in (2.13), and thus along the line a + (3 — 1 (and only this line) 
the mean- field theory is exact. 

2.3.2. Nonequilibrium partition functions and phase transitions In Sec. 3.3 we will 
show how the reduction process just described can be applied systematically to evaluate 
the matrix expressions (2.37) and (2.38) for the density and current in the ASEP. Here 
we quote the result for the normalisation Z^v 



from which the current follows via (2.42). We now discuss aspects of this expression 
which have generated recent interest. 

The relationship (2.42) between the current and the normalisation is intriguing. 
We may interpret Zjv as the partition function for the nonequilibrium steady state. 




(2.47) 




(2.48) 
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Then assuming that Zf^ grows exponentially with N we find that the analogue of the 
thermodynamic free energy is 

_ lim = hm In J . (2.49) 

Thus the current plays the role of the free energy. 

The phase transitions in the plane, manifested by different forms for the 
current J(q;,/?), arise from different large N asymptotic forms for Z^r and we shall 
determine these explicitly in Section 3. Moreover, recognising the normalisation 
as a nonequilibrium partition function, one can analyse the Lee- Yang zeros of the 
partition function and show that they predict precisely the location and order of the 
phase transitions in the driven system [17,61]. The interesting point is that the zeros of 
the partition function arc in the complex plane of transition rates as compared to 
equilibrium problems where partition functions zeros arc in the plane of some intensive 
thermodynamic variable such as fugacity or temperature. 

As we will discuss in Section 3.6, one can map the normalisation (2.48) of the ASEP 
to the partition function of a surface in thermal equilibrium with a heat bath. In this 
mapping, the transition rates a and /3 become equilibrium fugacities, and the phase 
transitions in the driven system correspond to adsorption-desorption transitions in the 
equilibrium interpretation. 

The expression for the normalisation (2.48) involves combinatorial coefficients 
known as ballot numbers [62]. There has been considerable recent interest in developing 
a combinatorial understanding of the origin of these numbers and other aspects of the 
steady state of the ASEP [63-66]. 

2.3.3. Comparison between matrix product states of equilibrium and nonequilibrium 
systems In the latter sections of this review we shall discuss in detail other 
nonequilibrium systems whose steady states are of matrix product form. For systems 
with open boundaries the ansatz used is of the form (2.35) where the number of matrices 
equals the number of possible states for each site. For systems with periodic boundary 
conditions a trace operation is used to obtain a scalar from the matrix product rather 
than using boundary vectors, see Section 4. At this point it is useful to establish what the 
difi^erences arc between these nonequilibrium steady states and the equilibrium steady 
states which may also be written in matrix product form. To illustrate this we rewrite 
the familiar transfer matrix solution of the the Ising model (an equilibrium system) as 
a matrix product state. 

Recall the one-dimensional Ising ferromagnet with N spins, Sj = ±1 and 
Hamiltonian 

= . (2.50) 

To this chain let us add two fixed spins, sl to the left of site 1 and sr to the right of site 
N . The Boltzmann weight of a configuration within this system when at equilibrium 
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with a heat bath at inverse temperature (3 can be written as a matrix product analogous 
to that used above for the ASEP (2.35). It takes the form 

f(si,S2,...,SN;SL,SR) =e~^^= {sl\Xs-,Xs2---XsJsr) (2.51) 



where 



and 



e' 



e-P' \ / 



(SL = +1| = ( 1 ) , (sL = -1| = ( 1 ) , (2.53) 

= = ( l-l^ ) ' = = ( 'e'/^"' ) ■ ^2.54) 
Note that the algebraic relations obeyed by the matrices, such as 

W{Si-i S[ Si+i Si-i Si Si+i)Xs,_iXs'.Xs,+i , (2.55) 

where W{Si_i Si Si^i — >• Si^iS^Si^i) is a transition rate for flipping spin Si, do not 
imply any reduction relations, rather they simply state the condition of detailed balance 
in the steady state. 

The partition function is given by 

Z = {sl\C''\sr) (2.56) 

where C = X+i + Calculating Z, for example by diagonalising C, one finds the 

free energy per spin 

Here we notice that the boundary conditions contribute only sub-extensively to the free 
energy, and that the free energy is an analytic function of temperature — that is, there is 
no ordering phase transition in the Id equilibrium Ising model. This contrasts strongly 
with the driven nonequilibrium steady state of the ASEP, where varying the boundary 
conditions induces phase transitions in bulk quantities such as the density. 

Mathematically this comes as a consequence of the D and E matrices for the ASEP 
being infinite-dimensional whilst the Id Ising transfer matrices are finite-dimensional 
and, because their elements are Boltzmann weights, non-negative. For matrices of this 
latter class one knows from the Perron-Probenious theorem that the largest eigenvalue 
is nondegenerate [67-69]. This means that is not possible, as temperature is varied, 
for the largest and second-largest eigenvalues to cross and the functional form of the 
dominant contribution to the partition function to change, which would give rise to 
a nonanalyticity in the free energy. This state of affairs applies, in fact, to all Id 
spin systems with short-range interactions [70]. On the other hand, if the interactions 
become long-ranged, or one goes into two or more dimensions, the transfer matrices 
become infinite-dimensional and phase transitions become a possibility. 
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To reiterate, the fact that one is with the matrix product construction solving a 
master equation, rather than simply packaging Boltzmann factors in an attractive way, 
means that one-dimensional nonequilibrium steady states can show much richer and 
more interesting behaviour than their equilibrium counterparts. As we shall see in the 
course of this article, the matrices that appear in various nonequilibrium models often 
violate the conditions for the Perron-Probenious theorem to hold. 



2.3.4- Comparison between matrix product states of quantum spin chains and stochastic 
systems Historically, matrix product states first appeared in the context of quantum 
spin chains [71,72]. There, the wavefunction |^) for the quantum spin chain is written 
as a product of matrices. Matrix product states also serve as a starting point for 
approximations and variational approaches such as Density Matrix RG [73,74]. Recently, 
matrix product states have received a lot of attention in the context of quantum 
information as they are able to describe highly entangled states [75]. 

The connection between the formalism for stochastic systems and that for quantum 
systems has been clearly reviewed in [76] . Here we summarise some important features. 
We first note that the master equation (1.1) can be written in a suggestive way as 

||P(t)) = -H\P{t)) (2.58) 

where \P{t)) is the vector of probabilities and if is a matrix formed from the transition 
rates. Clearly one can think of if as a Hamiltonian, \P{t)) as a wavefunction and (2.58) 
as a Schrodinger equation — this is referred to as the 'quantum Hamiltonian formalism' 
for stochastic systems (see e.g. [32,34]). The steady state of the stochastic system 
becomes the ground state of the quantum system and has eigenvalue 0. However, it is 
important to note some distinguishing features of (2.58) which are consequences of its 
true stochastic nature. First, the elements of \P{t)) are probabilities therefore must all 
be positive (this contrasts with the true quantum case where the amplitude squared of 
the elements are probabilities). Also, as the master equation conserves probability the 
sum of each column of the matrix H adds to zero whereas in the true quantum case 
there is no such restriction. 

For both a quantum spin chain and a nonequilibrium system the Hamiltonian can 
be written as a sum of local terms, for example representing the exchanges between 
neighbouring sites, 

ii = ^/i,,i+i. (2.59) 

i 

Now in the matrix product states constructed for quantum spin chains (see e.g. [77]) 
one finds 

hi,i+i\^) = (2.60) 

for each local operator For a stochastic system, the simple mechanism (2.60) 

to obtain a zero eigenvalue for H and hence the steady state, is only possible in the 
case of detailed balance, due to conservation of probability for each So to obtain 
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Figure 11. Dynamics of a system of first-class particles (circles) and second-class 
particles (pentagons). 

a true nonequilibrium steady state a more complicated mechanism is required. It is 
this more general mechanism which is at the heart of the matrix product formalism 
for nonequilibrium systems and gives rise to the algebraic relations, such as (2.39-2.41), 
which are a main focus of this review. 

2.4- Shocks and second-class particles 

So far we have considered the ASEP with only type of particle. In this section we discuss 
systems with two types of particles: first and second-class particles. The solution of the 
steady state of a system containing second-class particles was the second major success 
using the matrix product approach [78]. We shall discuss in detail the matrix product 
solution in Section 4, here we motivate the study of second-class particles by discussing 
some physical properties of interest. 

A system containing first-class particles (denoted 1) and second-class particles 
(denoted 2) has dynamics 

10-^01 20^02 12^21, (2.61) 

all processes occurring with rate 1 (see Figure 11). Thus a first-class particle treats the 
second-class particle as vacancy but from the point of view of a vacancy the second-class 
particle behaves like an ordinary (first-class) particle. 

Let us consider the velocity of a second-class particle in an environment of density 
p of first-class particles. The velocity of the second-class particle, denoted v, is given 
by the probability of the second-class particle hopping forward minus the probability of 
it hopping backward. Ignoring correlations between the second-class particle and local 
density the velocity is given by the density of a vacancy (in front), 1 — p, minus the 
density of first-class particles (behind), p 

v = l~2p. (2.62) 

Beware that this approximation is not at all correct! — there is structure in the local 
density profile around the second-class particle, as we shall see in Section 4. However, 
the expression (2.62) turns out to be exact and was proved by Ferrari [79]. We shall 
understand this result below. For the moment, note that (2.62) is precisely the velocity 
of a kinematic wave of density p discussed in Section 2.2. Thus the second-class particle 
travels along with the kinematic wave. 

As described in Subsection 2.2, when two kinematic waves of different densities 
meet a shock may form. The velocity is 1 — pi — p2 where pi is the density to the left 
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Figure 12. Illustration of coupling: the evolution of two systems one with an excess 
particle. The vertical lines indicate the bond which is updated in the two coupled 
systems. The dynamics of the excess particle (pentagon) is equivalent to that of a 
second-class particle — see text for explanation. 

and p2 is the density to the right of the shock and P2 > Pi- A second-class particle 
will be attracted to the position of the shock since when it is to the left of the shock 
1 — 2pi > 1 — pi — p2 and it catches up with the shock and when it is to the right of 
the shock 1 ~ pi — p2 > 1 — 2pi so that the shock catches up to it. Therefore we expect 
the second-class particle to track the shock. Moreover, one can use the second-class 
particle to actually define the position of the shock at the microscopic level. Then the 
density profile as seen from the moving frame of the second-class particle can describe 
the microscopic structure of the shock [80-82]. 

On an infinite system exhibiting a shock, it has been shown that the density as 
viewed from the second-class particle approaches a steady state distribution with density 
p+ at -|-oo and density p_ at — oo [79,81]. The exact structure of this steady distribution 
has been calculated using a matrix product approach on an infinite system which we 
will review in Sections 4 and 7.2.3. It demonstrates that the density profiles decay to 
their limits exponentially quickly in the distance from the second-class particle. An 
interesting case is when p+ = p_ so that there is a uniform density of particles and 
no apparent shock. However, in the moving frame of a second-class particle there is 
structure in the density profile: in front of the second-class particle the density is higher 
and decays to the asymptotic density according to a power law in the distance from the 
second-class particle; similarly behind the second-class particle the density is reduced 
and increases to the asymptotic density according to a power law. 

Another interesting aspect of the second-class particle is how its dynamics are 
related to the spreading of excess mass. The central idea, termed coupling [47] , is well- 
known in the mathematical community but less so within physics. Consider two systems 
containing only first-class particles, identical except that one system has M particles 
and the other M — 1 particles (one can consider either a finite system or an infinite 
system). The two systems start from initial conditions differing only by the position of 
the extra particle in the system with M particles. In order to implement the dynamics 
one can consider at each time step randomly choosing a pair of sites i, z + 1 to update; 
then if there is particle at site i and a hole at site i + 1 the particle is moved forward. 
In the dynamics let us choose the same pairs of sites in the two systems at each update 
(one can think of using the same random numbers in a Monte Carlo program). 

It is easy to convince oneself that after any length of time the configurations of 
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the two systems will differ only by the position of the extra particle (note that if we 
label the particles, the label of the extra particle will change under the dynamics). To 
see this consider the situation where the extra particle is at site i (see also Figure 12). 
Then, if any bond except i — or + l is updated, the result will be the same for the 
two systems and the extra particle will remain at site i. If site i + 1 is empty and bond 
+ l is updated then the extra particle at i will hop forward. If site i — 1 is occupied 
and bond i — 1, i is updated then the particle at i — 1 will hop forward in the system 
with M — 1 particle but not in the system with M particles. This results in the extra 
particle effectively hopping backwards. The latter two events result in the dynamics of 
the extra particle and the dynamics are precisely those of a second-class particle. 

Thus the dynamics, of a system comprising M — 1 first-class particles and one 
second-class particle describes the motion of an extra particle added to a system of 
M — 1 particles and the dynamics of second-class particles tells us about the dynamics 
of mass fluctuations. 

This idea can be used to recover the velocity (2.62) [78,81]. Consider adding an 
extra particle to a system of density p, increasing the density by Ap. The overall current 
increases from p(l — p) to p{l — p) + Ap(l — 2p) + 0{Ap'^). Now, since the system can 
be thought of as a system of first-class particles and an excess particle which behaves 
as a second-class particle we deduce that on a large system where Ap — > the speed of 
the second-class particle must be given by (2.62). 

An interesting result related to the spreading of mass fluctuations is that the second- 
class particle exhibits superdiff'usive behaviour. That is, with yt dcflncd as the distance 
travelled by the second-class particle, on an infinite system the variance in yt grows 
as [83] 



On a finite periodic system of size N, with a density p of first class particles, an exact 
calculation has shown [84] 



Thus on a finite system the behaviour is diffusive but with a diffusion constant which 
diverges with system size. 

These results are consistent with approximate calculations such as mode coupling 
[85, 86] for the spreading of excess mass fiuctuations which predict that the drift speed 
is 1 — 2p as the spreading of density fluctuations around the drift grows as t'^^^ on an 
infinite system. 

The dynamics of a second-class particle (2.61) is essentially passive in that it does 
not interfere with the dynamics of first-class particles. In this way the second-class 
particle can be though of as a 'passive scalar' for the TASEP. However, the density 
profile as seen from the second-class particle does have structure as explained earlier 
and gives insight into the structure of shocks. More recently generalisations of the 
second-class particle idea have been discussed in the context of passive scalars in other 
driven diffusive systems [87,88]. 



(2.63) 




(2.64) 



30 



1 a 1 

r~\ 






Figure 13. Dynamics of a system of first-class particles (circles) and defect particles 
(pentagons). 



2.4-1 ■ Defect particle The dynamics of the second-class particle can easily be 
generalised to hop rates 

10^01 20^02 12^21. (2.65) 

In the general case {3^1 the dynamics of a 2 particle is no longer passive and does 
affect the dynamics of the normal particles. We will refer to the dynamics (2.65) as 
that of a defect particle. A related defect particle dynamics has been considered as an 
exactly solvable model of two-way traffic problem involving cars and trucks [89]. 

An extreme case is /3 = where a normal particle cannot overtake the defect. In 
that case when a < 1 the defect particle is slower than the other particles and a 'traffic 
jam' may build up behind the defect particle. Actually there is an interesting phase 
transition which occurs on increasing the density: for low-density one has the 'traffic 
jam' phase where the density behind the defect particle is large and in front of the defect 
particle is small. For sufficiently high overall density there is a transition to a uniform 
density through the system. This jamming transition is mathematically analogous to 
Bose condensation, as we shall discuss in Section 8. 

For general a and (5 the steady state has been calculated using the matrix product 
approach [90]. A phase diagram emerges for a large system with density p comprising 
four phases (see Figure 22) characterised by different steady-state velocities, and density 
profiles surrounding the defect particle. We shall discuss this phase diagram in detail 
in Section 4.4. 



2. 5. The remainder of this review 

In this section we have introduced the asymmetric simple exclusion process (ASEP) and 
some of its relatives, and through various approximate treatments (kinematic waves, 
mean-field theory and the extremal current principle) have got a feel for some of the 
interesting physics, such as shocks and phase transitions, that might emerge. We have 
also presented the basic ideas underlying the exact matrix product approach that allows 
the existence of these phenomena to be confirmed. In the remainder of this review we 
fiesh out these ideas considerably, presenting detailed accounts of the various methods by 
which matrix product expressions can be evaluated and the physical properties of these 
systems fully revealed. As the review progresses, we will discuss all the main models 
that have so far been solved using the matrix product method. Mostly these models are 
multi-species generalisations of the ASEP in which particle number is conserved by the 
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Figure 14. Transitions into (dotted arrows) and out of (solid arrows) a sample 
configuration for (a) the ASEP on a ring and (b) the ASEP with open boundaries. 
Transitions occur at unit rate unless otherwise indicated. 

dynamics (except, possibly, at the boundary sites) whilst a small class of models further 
admits non-conservation in the bulk. Meanwhile, some of these models can be found by 
means of a systematic search under certain assumptions on the existence of reduction 
relations for the matrices similar to (2.39)-(2.41). Others involve more complicated 
structures which allow, for example, treatment of discrete-time updating schemes (such 
as parallel dynamics). Although in principle the matrix product method can be used to 
solve a very wide class of nonequilibrium stochastic dynamical models in one dimension, 
solutions of many interesting systems have proved elusive. We therefore conclude our 
work by highlighting a few of these models in the hope that, after reading this review, 
the reader will be inspired to tackle these new and challenging problems. 

For a full list of topics covered in sequence we refer the reader to the Table of 
Contents presented at the start of this review. We hope that this will also allow the 
reader to locate easily any sections that may be of particular interest. 

3. Detailed analysis of the ASEP with open boundaries 

Having introduced the ASEP with open boundaries (see Figure 3) and the matrix 
product expressions for the distribution (2.35), density (2.37) and current (2.38) in the 
steady state we now explain in detail how these are calculated using the matrix reduction 
relations (2.39)-(2.41). First, though, we must demonstrate that the distribution 
implied by (2.35) and (2.39)-(2.41) is, in fact, the stationary solution of the master 
equation for the process. 

3.1. Proof of the reduction relations 

3.1.1. Domain-based proof As a first step towards understanding why the matrix 
product solution of the ASEP with open boundaries given above is correct, we consider 
the simpler case of the ASEP on a ring. A typical configuration is shown in Figure 14a. 
Recall that the master equation (1.1) for the probability of being in configuration C has 
two contributions: a gain term from transitions into C and a loss term from transitions 
into another state. In the steady state, we require these loss and gain terms to balance. 
As is evident from Figure 14a there is for each domain of particles on the ring one 
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way to enter a given configuration C (each has a particle joining the end of a domain) 
and one way to exit (when a particle leaves the front of a domain). Therefore we have 
for any given configuration C 

Y^W{C' ^C)^Y.W{C^C') . (3.1) 

c c 

The loss and gain terms in the master equation (1.1) thus cancel if P{C) = const VC. 
Putting this another way, every gain term P{C')W{C' C) exactly balances one of the 
loss terms P{C)W{C — > C"). Note this is similar, but not the same, as the cancellation 
that occurs when detailed balance (1.2) is satisfied since C ^ C". This more general 
cancellation scheme is sometimes referred to as dynamic reversibility [9] or pairwise 
balance [91]. 

When one has open boundaries, as in Figure 3, the number of ways into a given 
configuration does not always equal the number of ways out. For example, when the 
left boundary site is occupied and the right boundary empty, the number of ways in 
exceeds the numbers out, and furthermore not all hops occur at the same rate (see 
Figure 14b). The strategy to finding the steady state solution of the ASEP is to have a 
partial cancellation between the terms corresponding to particles attaching to the rear 
or detaching from the front of domains in the bulk. The remainder is then cancelled 
with terms coming from the boundary interactions. As we now describe by reference to 
particular examples, it is this cancellation scheme that is implied by the matrix relations 
(2.39)-(2.41). 

For our illustrative purposes, let us consider a particular six-site configuration 
_Q0_0_ for which the master equation reads 

^/(_Q0_0_) = - a/(lQO_0_) 



+ /(61o_o_)-/(_o61o_) 



+ /(_Qo61_) -/(_oo_6:' 
+/9/(_Qo_oo) . (3.2) 

In the gain terms (i.e., those with a positive sign), the arrows indicate the hop that 
leads to the configuration on the left-hand side; in the loss terms, the arrows indicate 
the transitions by which this configuration can be exited. On the right-hand side we 
insert now the matrix product expressions: 
d 

^/(_QQ_Q_)= -a {W\E DDEDE\V) 

+ {W\DEDEDE\V) - {W\EDDEDE\V) 

+ {W\EDDDEE\V) - {W\EDDEDE\V) 

+ P{W\EDDED D\V) . (3.3) 

Here we have underlined the mathematical objects (vectors or matrices) that correspond 
to the sites or boundary reservoirs involved the transitions indicated with arrows in (3.2). 
We note that all the these underhned terms are of the form (W\E, DE or D\V) and so 
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can be reduced using (2.39)-(2.41). Using these reduction relations we find 

^/(_Q0_0_)= - {W\DDEDE\V) 

+ {W\DDEDE\V) + {W\EDEDE\V) - {W\EDEDE\V) - {W\EDDDE\V) 
+ {W\EDDDE\V) + {W\EDDEE\V) - {W\EDDEE\V) - {W\EDDED\V) 
+ {W\EDDED\V) , (3.4) 

where here the underhned vectors and matrices show what remain after using the 
reduction relations. We see that in the second and third lines of this expressions, 
which originated from particles joining and leaving domains, the middle two terms 
cancel; furthermore, in both cases the first term cancels with the last term on the 
previous line, and the last term with the first term on the next. Hence the weight 
/(_QO_o_) = (W\EDDEDE\V) is stationary, as claimed. If one has a configuration 
comprising n domains of particles, none of which contact either boundary, one finds a 
master equation of the form (3.2) with n pairs of terms like those in the second and 
third fines of (3.2), and the same cancellation mechanism goes through. For the case of 
the empty lattice, n = 0, one can verify that the term arising from the left boundary 
cancels with that from the right. 

It only remains to check that weights corresponding to configurations with particles 
on the left, right or both boundary sites are also stationary. Let us consider a 
configuration that has a domain of particles at the left boundary. For such a 
configuration, the master equation begins with the terms 

^/(QOO_ . . .) = «/(1qo_ . . .) - /(QQOI ...) + ... (3.5) 

= a {W\E DDE-- ■ \V) - {W\DDDE ■ • • |y) + • • • (3.6) 

= {W\DDE ■■■]¥)- {W\DDE- • • [F) - {W\DDD-- + . (3.7) 

The first two terms on the last line cancel, leaving a single negative term. One can 
also verify that this term will cancel one coming from the next domain of particles on 
the lattice, or from the right boundary. Similarly, if one has a domain of particles in 
contact with the right boundary, after using the reduction relations and performing 
a cancellation, a single positive term remains that cancels with one coming from the 
previous domain of particles on the lattice, or from the left boundary. Finally, after 
checking that the matrix-product expression for the statistical weight of fully occupied 
lattice is also stationary, we can conclude that the matrix product weights, defined 
by Equation (2.35), in conjunction with the reduction relations (2.39)-(2.41), give the 
steady state solution of the ASEP, as previously claimed. 

3.1.2. Algebraic Proof We now show how this pictorial representation of the 
cancellation mechanism between domains can be encoded in an algebraic way. This 
will pave the way to understanding the general results of [92] that we review in Sec. 9. 
The master equation can be written as a sum of terms, each relating either to a 
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pair of neighbouring sites, or to one of the end sites, i.e., as 



dt 



f{n,T2,. . . ,Tn) 



N-1 



i=l 



/(ri,r2, ...,rjv) 



(3.8) 



The operators h construct the required gain and loss terms in the master equation via 

hi,i+lf{- •■ ,Ti, Ti+i, • • •) = W{Ti+iTi TiTi+i)f{- ■ ■ , Ti+1, Ti, ■ ■ •) 

-W{TiTi+i Ti+iTi)f{- ■■,Ti, Tj+i, ■ ■ ■) (3.9) 

hLfiri, • • •) = Wi{l - Ti ^ Ti)/(1 - Ti, • • •) - Wi{ti ^ 1 - Ti)/(Ti, • • •) (3.10) 

flRfi- • • , Tiv) = Wn{1 -Tn^ TN)f{- • • , 1 - Tiv) - W^iv(TAr ^ 1 - Tiv)/(- ■ ■ , Tn) (3.11) 

in which W{tt' — > r'r) gives the rate at which a neighbouring pair of sites in 
configuration (r, r') exchange, Wi{ti — > 1 — Ti) gives the rate at which site 1 changes 
state and Wn{tn — > 1 — tjv) gives the rate at which site N changes state. For the ASEP 
we have 

iy(10^01) = l (3.12) 
Wi(0^1) (3.13) 
Wn(1^0) =/3 (3.14) 

and all other rates zero. 

In terms of matrix product expressions (2.35) for the weights, (3.9-3.11) become 

k,+i{W\ ■ ■ -X^^^AXnX^.jX,.^, ■■■\V) = {W\ ■ ■■X,^_,[W{n+in ^ r,r,+0^r,+i^n 

- Winn+i r,+iri)X,^X,^jX,^_^, ■ ■ ■ |\/) (3.15) 
hL{W\Xr, ■■■\V) = {W\ [Wi{l - n ^ ri)Xi^r, - Wiin ^ 1 - ri)Xr,] ■■■\V) (3.16) 
hR{W\ ■ ■ ■ [Wn{1 -tn^ rN)Xi.r^ - W^irN ^ 1 - r^)X.^] \V) (3.17) 
where we have stretched the notation so that, for example, now acts on the matrices 

Xt- , X-r- I -, . 

Meanwhile, we suppose there exists a set of auxiliary matrices Xi such that 
additionally 



k,+i{W\---X,^X^^^ 
hL{W\X,^---\V)^ 



■■\V) = {W\---[Xr^X 

-{w\xJ...\v) 



Ti+1 



XT^Xr^_^_l] 



hR{W\---XrjV) = {W\--- XrjV) 



\V) (3.18) 
(3.19) 

(3.20) 



If this is the case then the sum in (3.8) telescopes to zero (i.e., there is a pairwise 
cancellation of terms). Therefore, if the right-hand sides of (3.15)-(3.17) and (3.18)- 
(3.20) can be shown to be equal, it follows that a stationary solution of the master 
equation for the process has been determined. Equating the terms in square brackets 
in (3.15)-(3.17) and (3.18)-(3.20) is sufficient to ensure overall equality of the two 
expressions. That is, if 



W(Ti+iTi TiTi+i)Xri^^Xri - W (jiTi+i Ti+iri)X^.X^^_^^ = X^.X^^,^^ - Xr^Xr^^^ (3.21) 
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w,{i - n ^ n){w\x,^r, - w,(n ^ i - n){w\x,^ = -{w\Xr, (3.22) 

Wn(1 - tat ^ TAr)Xi_,^|y) - WNiTN ^ 1 - rN)XrJV) = . (3.23) 

Note that (3.21-3.23) amount to a total of eight conditions which become with the 
set of rates (3.12)-(3.14) 

^ DD - DD ^ EE - EE (3.24) 

DE =ED-ED = -DE + DE (3.25) 

a{W\E = - {W\b = {W\E (3.26) 

pD\V) = E\V) = -D\V) . (3.27) 



With a scalar choices for the auxiharies, viz D = —1, E = 1 and ^4 = 0, these reduce to 
the three equations (2.39-2.41) as required. 

The algebraic proof can be easily generalised to other models, for example those 
with more general bulk dynamics than exchanges, as shall be carried out in Section 9. 

3.2. Calculation of the nonequilibrium partition function 

Our first calculational task is to compute Zn given by (2.36), from which the current 
follows via (2.42). There are several approaches to the task and these are representative 
of the three main approaches to calculation which will be employed at different points 
in this review: direct matrix reordering, diagonalisation of an explicit representation of 
the matrices and generating function techniques. 

3.2.1. Direct calculation by matrix reordering In Section 2.3.1 we explained how 
the reduction relations can, in principle, be used repeatedly to bring any arbitrary 
matrix product U into a standard form (2.45) which can then be straightforwardly 
evaluated (2.46). We now put this principle into practice, and begin with the product 
= {D + E)^ required to calculate the normalisation (2.36) and thence the current 
(2.38). By performing the reduction longhand for the first few N, one finds 

C =D + E, (3.28) 
C'^ ^ DD + DE + ED + EE = DD + ED + EE + D + E , (3.29) 
= DDD + DDE + DED + DEE + EDD + EDE + EED + EEE 
= DDD + EDD + EED + EEE + 

DD + DE + DD + ED + DE + EE + ED + EE 
= DDD + EDD + EED + EEE + 2{DD + ED + EE) + 2{D + E) . 

(3.30) 

Eventually one is drawn to the general formula 

«" = E^f^^t^'^-'- (3.31) 
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which can be proved by induction [58]. For this two prehminary identities are required. 
First, 

D^C = D + + • • • + + E , (3.32) 

which can be proved inductively by multiplying both sides from the left by D and using 
(2.39). The second identity involves the combinatorial factors 



B 



N,p 



N](N-p)\ 







otherwise 



It reads 



Bn,p — Bn,p+i 



(3.33) 



(3.34) 



V7V>l,p>0 

which can be verified by substituting in the explicit expressions from (3.33). 

We now assume that the identity (3.31) for has been shown to hold for some N, 
and derive an expression for (7^+^ by multiplying from the right with C. Then, (3.32) 
implies 

TV p 

^ J2 Bn^p [^'^' + E'^iD + + ■ ■ ■ + Df-'^+i)] . (3.35) 

p=0 q=0 

Now, using (3.34) we find 

N+l p-1 

= ^ B^^^^p [E'i+' + Ei{D + D' + --- + DP-'i)] - 

p=l q=0 
N+l p-2 

J2 Bn+i,p J2 [^'^' + ^'(^ + D^ + --- + DP-'^-')] (3.36) 

p=2 q=0 

n+l 

= J2 ^N+i,p [EP + EP-^D + EP-^D^ + --- + DP] (3.37) 

p=i 

which is equal to the right-hand side of (3.31) with N —>■ N + 1. Since (3.31) is obviously 
correct for = 1 (C = D + £■) , it follows by induction that it is true for all A'^. Using 
(2.40), (2.41) and the convention = 1 we immediately find the exact expression 



N 



p=Q 



q=0 



a 



N 



p=0 



i _ i 

a 



(3.38) 



3.2.2. Diagonalisation of an explicit representation So far we have not given any 
exphcit representations of the matrices D,E and vectors |y) that satisfy (2.39- 
2.41). Several possible representations were proposed in [58]. For example, one can 
take: 

/1100---\ /1000---\ 



D = 



110 
11 
1 



V : 



E 



110 
110 
11 



(3.39) 
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{W\ = k(^1, a, a? . . ) \V) = k 



( 1 \ 

h 

V • I 



(3.40) 



where 



1 — a 



a — 



6 = 



a. 



1-/3 
/3 



(3.41) 



= (a + /3 — 1)/q;/? is chosen to ensure that V) = 1. Note that for this we require 
q; + /3 > 1, but as we shall discuss later (in Section 5.3) one can analytically continue 
results to q; + /3 < 1 in a straightforward manner. Some other representations are 
presented in Appendix D. 

Let us consider the eigenvectors of C, given in this representation by 

/ 2 1 . . \ 

12 10 



C 



12 1 
12 



(3.42) 



V • • / 

We parametrise the eigenvalues by 2(1 + cos 6*) where < 6* < 27r and denote the 
associated eigenvector by | cos ff) so that 

C|cos^) = 2(1 + cos^)|cos^) (3.43) 

where 

/ sin^ \ 
sin 26' 

sin 3^ . (3.44) 



cos( 



sin^ 



V 



/ 



Note that the elements of | cos Q) arc precisely the Chebyshev polynomials of the second 
kind [93]. It is easy to show (by summing geometric series) that 



(iy|cos^) = ^ J]a"sin(n + 1)^ 



n=0 



(1 - ae^^)(l - a^-'^) 



(3.45) 



and (cos 6^1 y) is given by a similar expression with a replaced by h. Also it is easy to 
show (using orthogonality of sin nQ for different n) that 

l = i/ desin2^|cose)(cos^| . (3.46) 
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Therefore Z^^ can be evaluated as follows 

Zjv ^l^\C^\V) (3.47) 
1 /■^'^ 

= -/ d^sin2^(PF|C^|cos^)(cos^|V) (3.48) 

= ^ sin^g[2(l + cosg)]^ 

TT 7o (1 - ae^^)(l - ae-^«)(l - &e^^)(l - 6e-^^) ' ^ ^ 

This integral representation of the normalisation can be shown to be equivalent to (3.38) 
by using the calculus of residues (see Appendix C). 

3.2.3. Generating junction approach Let us define the generating function of the 
normalisation as 

oo 

Z{z) = ^""^^ ■ (3-50) 

In order to determine this quantity in a fast way the idea is to work with formal power 
series involving matrices, i.e., expressions of the form 

oo 

/([/) = J] a„C/" (3.51) 

n=0 

where C/ is a matrix. The "formal" component of this approach lies in the fact that 
we do not worry about convergence properties of these scries. Rather, we simply use 
such expressions as a device for manipulating matrix product expressions. For example, 
one might multiply both sides by some other formal series in a matrix V ^ bearing in 
mind that if U and V are noncommuting, one must multiply both sides from the same 
direction. 

Consider the formal series 

oo 

= y z"C" (3.52) 



l-zC ^ 

n=0 

SO that the generating function of the normalisation (3.50) can be written as 

2{z) = {W\y^\V) . (3.53) 
As was observed by Depken [94], Equation (2.39) has the simple consequence that 

(1 - r]D){l -r]E) = l- r]{D + E) + rj^DE = I - r]{l - r])C . (3.54) 

Putting z = ri(l ~ rj), and multiplying both sides from the left by i-^j from 
the right by learns that 

1 _ 1 1 
1 - zC ~ l-r]El-r]D ' 
Using (2.40) and (2.41) one then finds that 

Z(.)=fl-!^V7l-!^) . (3.56) 



(3.55) 



a J \ (3 
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Note that the condition Z{0) — 1 imphes 77(0) = which means that we must choose 
the root 

rj{z) = ^ (1 - VT^) . (3.57) 

The generating function approach provides a very fast route to the phase diagram 
and current via the asymptotic analysis of (3.56) that we perform Sec. 3.4. Furthermore, 
the exact expression (3.38) can also be recovered with very little work through an 
application of the Lagrange inversion formula for generating functions which we recap 
in Appendix B. As indicated in Appendix B, this formula can be applied because the 
combination ri{z) / z can be expressed in terms of ri{z) alone, i.e., its dependence on 
z enters only implicitly through the argument of rj^z). Then, the Lagrange inversion 
formula implies that the normalisation Z]\f is given by 



1 / <9 




N 



= -^l"ai + ''sifj(""'l3|i3|(i3^j (3.89) 

where we introduce an important notation: 



{x'^}f{x) denotes the coefficient of x" in the formal power series f{x) 



Expression (3.60) is in perfect agreement with (3.38), as required. 

3.3. Exact density profile from direct matrix reordering 

In order to calculate the density profile we use the approach of matrix reordering. This 
has the advantage of giving an exact finite sum expression (3.65) for finite systems. We 
start from (2.37), which contains the matrix product DC^ . Again, one can perform the 
reduction manually for small j (we shall omit this here) and deduce the identity 

i-i i+i 
DC^ = J2 Bk+i^iC'-" + Yl B,,k-iD^ . (3.61) 

fe=0 k=2 

This was proved in [58], again using an inductive argument. Assuming (3.61) holds for 
some J > 1, we find this implies for j ' + 1 that 

i-i i+i 
{Da)C = J2 Bk+i,iC^^'-' + Yl Bj,k-iD'C (3.62) 

k=0 k=2 

i-1 i+i 
= J2 Bk+i,iC^+'-'' + ^(S,+i,fc - Bj+,,k+i){C + D^ + D' + --- + £''=+^3.63) 

k=0 k=2 
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Figure 15. Particle current in the ASEP at finite system size for two different 
combinations of a and (3. 

where both (3.32) and (3.34) have been used. Since, for > 1, i?Ar+i,i = a 
httle more algebra reveals that 

DC^+i = ^ + + 5,+i,2D3 + ■ ■ ■ + (3.64) 

which is (3.61) with j ^ j ' + 1 as we desired. The case j = 1 is easily verified using the 
reduction relation (2.39), which completes the proof of (3.61) for all j. An expression 
for the density then follows by substituting (3.61) into (2.37). We find 

N-i ^ ^ N-i ^ 

= T.Bn,^ + i^E^^--^ (3-65) 

n=l p=l 

for 1 < i < A^. Note that the case i = A^ is special: pn = Z^-i/ {[3Zm). 

In the next subsection, we shall analyse the large- A^ behaviour of the current and 
density profiles. First, though, we remark that the formulae (3.38) and (3.65) are of 
utility if one wishes to study finite systems, e.g., in the context of the biophysical 
applications of the ASEP where the thermodynamic limit is not appropriate. For 
example, the particle current as a function of system size is plotted in Figure 15 for 
the combinations a = [3 = 1 and a = l,/3 = |. The figure is suggestive that the decay 
to the asymptotic current with increasing system size is slower for the case a = (3 = 1 
compared to a = 1 , /? = | . As one might guess from the mean-field treatment of the 
ASEP in Section 2, these different decay forms arise from being in different phases of 
the model. 

3.4- Asymptotic analysis of the current and density profiles 

In order to obtain the exact phase diagram for the ASEP, it is necessary to determine 
the asymptotics (large- A^ forms) of the current and density profiles. There are a number 
of ways to achieve this. A method that is systematic and easy to apply in a wide range of 
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cases is to construct the generating functions functions of the desired quantities. Then 
one can apply standard techniques for extracting the asymptotics of the coefficients 
from a generating function (see, e.g., [95]). For convenience, we have recapitulated 
these techniques in Appendix B. 

It will also be useful to have at our disposal some generating functions for common 
combinatorial quantities. We first note that ri{z) (3.57) that appears in (3.56) is actually 
(up to a factor of z) the generating function for Catalan numbers: 

N=l ^ ^ 

The Catalan numbers C„ = ;^(^^) appear in a remarkably large number of 
combinatorial contexts [96] and will also be used many times in the course of this 
article. Equation (3.66) can be seen by expanding the square root in (3.57) in powers 

of 2;. It is also instructive to consider the generating function for ballot numbers 

00 

G^{z) = J2 (3.67) 

N=p 

where the coefficients B^^p, are defined in (3.33), and we have that Go{z) = 1 and 
^1(2;) = ri{z). It turns out that Gp{z) is simply the p^^ power of rj{z). To see this, note 
that (3.34) implies 

Gp{z) = Gp+i{z) + zGp.i{z) (3.68) 

and so that if Gp{z) — [r]{z)]P one must have r]{z)[l — r]{z)] — z which is indeed satisfied 
by (3.57). 

As noted in Appendix B, the asymptotics of the coefficients of a generating function 
(here, the normalisation Zn) are dominated by the singularity nearest the origin. Now, 
the function r]{z) has a square-root singularity at z — j which is always present. 
Additionally there is a pole at 77(2;) = a; however, this is present on the positive 
branch of the square- root that appears in 77(2;) only if a < |. Similarly, there is a 
pole at r){z) — (3 ii (3 < ^. As r){z) increases monotonically from zero in the range 
< 2; < |, the dominant singularity is the square- root singularity at 2; = | if both 
a and f3 exceed |. Otherwise, whichever of « or /5 is smaller supplies the dominant 
singularity — see Figure 16. Note that the regions within which particular singularities 
dominate coincide with the phases previously established using a mean-field approach, 
shown in Figure 5. In other words, this analysis of the generating function reveals the 
phase diagram to be exact. 

In the region identified as the low-density phase from the mean-field calculation, 
a < P,a < ^, the pole in (3.56) at zq — Oi{l — a) implies — via the results presented in 
Appendix B — that 



Zn{ol,(3) 



lim ( 1 - — ) Z{z-a,(3) 

z^zo \ Zq 



a(3 
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(3.69) 
(3.70) 



(i) 



(ii) 



(Hi) 



a(l-a) 



P(l-P) 



Figure 16. Poles (crosses) and branch cuts (solid lines) in the complex-z plane of the 
generating function Z{z;a,f3) at three different values of a and (3. (i) a > i, /3 > i; 
(ii) a < J, (3 > ^] and (iii) a < i, < i and (3 < a. In each case the singularity 
dominating the asymptotics of the normalisation is the one closest to the origin and 
indicated with an arrow. 



77, ^"(1-") • (3.71) 

Using (2.42), one finds the current is J = zq = a(l — a). In the High-density phase, 
/? < a, /3 < I the pole aX zq = (3{1 — [3) dominates and one obtains the same expression 
for the normahsation but with a and (3 exchanged. This is in accordance with the 
particle-hole symmetry of the model. In the maximal-current phase, where both a and 
l3 are greater than i, we expand up the square- root in (3.56) to find 

Ziz] a, = ( 1 - \^— + —7^1 v^r^4i + 0(1 - Az)] . (3.72) 

^' (2a- l)(2/5- 1) V [2a-l 2/3-lJ ^ J 

It then follows from the results of Appendix B (Equation B.7) that 

4a(3{a + f3-l) 4^ 
^^^"'^^^ 0F(2a- 1)2(2/3 -1)^A^ ^^-^^^ 
which can be shown to agree with the expression given in [58] after a little algebra. In 
this phase, the current J = zq = j, thus showing that the mean- field predictions for the 
asymptotic current in all three phases, as presented in Table 1, are exact. 

We now turn to the exact formula for the density profile (3.65) which for convenience 
we give again here: 

N-i N-i 
n=l p=l 

Although formidable at a first glance, its analysis is in fact reasonably straightforward 
when both N and i are large. More precisely, we shall take N ^ 00 and i ^ 00 with 
the distance j = N — i from the right boundary fixed. It is easily verified that in this 
limit one has in all phases 

lim = (3.75) 



N^oo 
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for any fixed finite j. Recall that J is the current the thermodynamic limit, as 
determined above. Thus 

toj^p._,=j:B,,j"+:^5:^. (3.76) 

The second sum in this expression we have met before: it is simply the normalisation 
for the ASEP (2.48) in the hmit a — > oo. Hence, for j » 1 we can use the asymptotic 
analyses from above to determine that 

(3{l - 2(3) 



3 R 



p=\ 



4' 



/34^' 



/3<i 

■ (3.77) 

The case (3 — special because the pole and square-root singularities in the generating 
function coincide, thus modifying the nature of the dominant singularity. 

The first sum can be analysed through its generating function in a manner similar 
to that described for the normalisation above. One ultimately finds 

j:B„,J^ = {y^}^ (3.78) 

where we recall that {y^} picks out the coefficient of from the generating function. 
In the high- and low-density phases, where J < |, the pole provides the dominant 
contribution to the sum for j' — > oo, with finite-j corrections supplied by the square- 
root singularity in r){y J) at y — jj. The result is 

J-^ T (A TV 

EBn..J'~v(J)-^j^ + -. (3.79) 

Note that r]{J) = mm{a, f3} in these phases. Meanwhile, in the maximal-current phase, 
the pole and square-root singularity coincide and then 

^ B^,,J^ ^Ul-^ + 0(i-^/^)) . (3.80) 

n=l ^ V VTTJ / 

3. 5. Exact phase diagram 

We can now put these results together to determine the density profile at a fixed (but 
large) distance j from the right boundary in the limit N ^ oo using (3.65). We do not 
need to consider the left boundary explicitly because of a particle-hole symmetry in the 
model that implies 

PN+i-i{a, P) = 1- pi{(3, a) (3.81) 

since holes enter the system at the right boundary at a rate f3, hop with unit rate to the 
left and exit at rate a. There are a number of different cases to consider in obtaining the 
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density profile, since not only does the first term in (3.74), and prefactor of the second, 
depend on which phase the system is in, but the second sum additionally takes two 
different functional forms at the right boundary in the low-density phase. The various 
distinct cases are as follows. 

• Maximal-current phase /9 > |,q; > |. Here, the generating- function analysis 
of (3.65) involves only algebraic singularities, rather than poles. Hence one finds a 
power-law decay of the the density from the bulk value of |, viz, 

1 ' 
2 



-^ + o(r^/2) 



(3.82) 



Although the mean-field theory predicts the bulk density correctly, the decay 
exponent is incorrectly predicted as 1 rather than i. Note also the universal (i.e., 
independent of a and /?) prefactor of this decay which is also found at the left 
boundary as a consequence of the particle-hole symmetry. 

Low-density phases a < < (3. As noted above, the second term in (3.65) 
takes different forms depending on whether /3 is less than or greater than |. This 
gives rise to two sub-phases. 

- LD-I /5 < |. Here, the generating function for the density contains only poles, 
and so a purely exponential decay is obtained at the right boundary: 



p^., ^ a + (1 - 2/?) [J^^^ ) ■ (3.83) 

— LD-II /3 > |. In this instance a pole and algebraic singularity combine to give 
an exponential decay modulated by a power law: 

~ «+ - (2^] v5^* '°<^ - ■(^•«^' 

High-density phases /3 < < a. In the high-density phases, all variation in 
the density is at the left boundary (which is out at infinity in the foregoing analysis) . 
At the right boundary the density takes the constant value 

~ 1 - /5 , (3.85) 
as was seen in the mean-field analysis. 

High/low-density to maximal-current transition lines Along the boundaries 
between the high- or low-density and the maximal- current phases, the bulk density 
is I with deviations of order (the prefactor can be calculated by developing 

the asymptotic expansions given above: this is left as an exercise for the reader). 
High- to low-density transition line Here, the situation is more interesting. 
Since a — f3 < ^, the generating function for the normalisation has a double pole, 
and so ~ (const)A'"J~^. This factor of N in the normalisation gives rise to a 
linear density profile from the prefactor Zi^xjZ^ of the second sum in (3.65). That 
is, 

A = a + (l-2a)^. (3.86) 
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Figure 17. Exact phase diagram for the ASEP. The two high- and low-density sub- 
phases differ in the decay of the density profile to its bulk value from the boundaries. 

This density profile can be understood from the study of shock fronts discussed 
in Section 2.4. A shock, at x = 0, is stationary if the density approaches Pl < \ 
as the spatial coordinate x — > — oo and a value p/j = 1 — pz, > as a; oo. 
One notes that this is the situation imposed by the boundary conditions on the 
open system in the N oo limit along this transition line: pl = a, Pr = 1 — a- 
One has a diffusing shock front separating a low-density region to the left and a 
high-density region to the right. Given that this shock front performs a random 
walk, its position on the lattice is given by a fiat distribution which in turn implies 
a linear average density profile of the form given above. 

With these results we are now able to present the complete phase diagram for 
the ASEP. It is shown in Figure 17 and differs from the mean-field version (Figure 5) 
in that the high- and low-density phases have acquired two sub-phases in which the 
density decay near the boundary takes a different form. The lengthscale characterising 
this density decay diverges as the maximal- current phase is approached from either the 
high- or low-density phases. This type of physics is suggestive of a continuous phase 
transition [59]. Meanwhile, the phase coexistence at the high- and low-density phase 
boundary is reminiscent of phase coexistence at a first-order phase transition in an 
equilibrium system. Motivated by these observations, we now explore more deeply the 
relationship with equilibrium phase transitions. 

3.6. Connection to equilibrium phase transitions 

In Section 3.4 we saw that the generating function of Catalan numbers played a key part 
in the analysis of the ASEP's phase behaviour. This fact can be used to make contact 
with equilibrium phase transitions in two-dimensional surface models [97-99] since one 
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Figure 18. A Dyck path is a random walk on the rotated square lattice that returns 

to the origin for the first time at a fixed end-point. Between the initial upward and 
final downward step, the walk comprises a some number (possibly zero) of Dyck paths 
placed end-to-end. 

of the many combinatorial objects counted by Catalan numbers are excursions on the 
rotated square lattice, known also as Dyck paths (see e.g., [100] for further applications 
and references). 

Specifically, a Dyck path of length 2N comprises N up-steps and N down-steps, so 
that it contacts the origin both at the start and the end. It is further constrained never 
to contact the origin in the intermediate steps — see Figure 18. The coefficient of in 
the generating function G{z) counts the number of paths of length 2N . Referring to 
Figure 18 we see that a Dyck path of a given length 2N begins with an up-step, ends on 
a down-step and comprises any number of excursions that together have length 2(A^ — 1) 
as can be seen from Figure 18. Therefore the generating function for Dyck path obeys 

G{z) = z[l + G{z) + G\z) + G\z) + ■■■] (3.87) 

where on the rhs of (3.87) the factor z comes from the pair of initial and final up and 
down steps and the sum of terms G'^{z) come from excursions which return exactly n 
times to site 1 in between. We see then from (3.88) that G{z) — rj{z), the function 
defined in (3.57). 

Similar reasoning implies that the the coefficient of z^ a~'^(3~'^ in the generating 
function (3.56) 

.W^(l-M£))-'(l-f)- (3.8.) 

counts the number of walks of length 2N that return to the origin a total of n -|- m times 
(the last of those returns being the end of the path). It is convenient to have the first 
n of these excursions take place above the origin, and the remaining m below. Then, 
the coefficient of z^ , which is the normalisation of the ASEP Z]^ for an A^-site system, 
is also the partition function for the equilibrium ensemble of surfaces of length 2N that 
is constrained to start and end at the same height, crosses the origin at most once (and 
then only from above), and has a fugacity 1/a associated with each contact with the 
origin from above, and with those from below — see Figure 19. This construction 
has been dubbed a one-transit walk [98]. 

It is an elementary exercise in statistical mechanics to derive the free energy per 
site in the thermodynamic limit as a function of the density of contacts pa from above 
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Figure 19. A one-transit walk whose statistical weight in the fixed-length ensemble 
depends on the density of contacts from above and below. 
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Figure 20. Phase diagram of the one-transit walk in the plane of equilibrium fugacities 
a and (3. 

and pb from below, these being the order parameters of the surface. This free energy 
is [98,101] 

f{pa, pb) = {pa \n.a + pb In P} 

-{{2 -Pa- Pb) ln(2 -Pa- Pb) -{1-pa- Pb) ln(l -Pa- PbW.90) 

in which the first term in curly brackets is an energy-like quantity that is lowered by the 
surface contacting the origin and the second an entropy-like quantity that is increased 
by the surface making large excursions away from the origin. The equilibrium state 
arising from this energy-entropy competition is found in the conventional manner, i.e., 
by minimising the free energy with respect to the order parameters pa and pb- 

The free energy is minimised by choosing pa, Pb to lie somewhere on the boundary of 
the physical region < pa + Pfe < 1- When a and /3 are large, the entropy of walks that 
rarely contact the origin dominates and the surface is desorbed: Pa = P6 = 0. On the 
other hand, if one of a or /5 is smaller than the critical value of |, contacting the origin 
is energetically favourable. If a < (3, it is most favourable for all contacts to come from 
above, otherwise from below. Thus the low-density, high-density and maximal- current 
phases of the ASEP correspond to two adsorbed phases, one from above and one from 
below, and a desorbed phase respectively — see Figure 20. We remark that a similar 
unbinding transition was explored using the mapping from ASEP to the disordered 
equilibrium problem of directed polymers in a random medium [102] . 
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An interesting feature of this association of nonequilibrium and equilibrium phase 
behaviour is that the orders of the phase transitions must agree in both cases. In the 
equihbrium case, these are defined in the conventional way, i.e., through discontinuities 
in derivatives of the equilibrium free energy with respect to the contact fugacities 
1/a, In the nonequilibrium case, on the other hand, derivatives of the "free energy" 
(the logarithm of the normahsation — see Equation 2.49) with respect to the transition 
rates a and (5 have no obvious connection to moments of physical observables as they 
do for an equilibrium system. We remark that this phenomenon had previously been 
observed in a study of the Lee- Yang zeros of the ASEP normahsation [17] and direct 
the interested reader to two recent thorough reviews of this subject [61, 103] for details. 



3.7. Joint current- density Distribution 

Finally, in this section on the ASEP with open boundaries, we develop the generating 
function approach further to compute certain more complicated quantities. It is useful 
to introduce a pair of matrices 

E'^(l-x- y)Y^ ■ (3.92) 

Then, for example, the combination ^_^_y D' in a generating function constructs a 
domain of particles, the size of which is conjugate to the parameter (fugacity) x. These 
matrices have the special property that they satisfy the familiar relation 

D'E' ^D' + E' (3.93) 

as a consequence of (2.39). This means that results already obtained by, for example, 
direct ordering of D and E matrices can be carried across to expressions involving D' 
and E'. One quantity calculated in this way is a joint current-density distribution was 
calculated in [104], the derivation of which we now briefly review. 
The interest here is in the quantity 

P^^M,K) . ^{.M.-MKH'I^ (-^-l^)-^ _J_|V.) . (3.94, 

This gives the probability that a configuration of N sites containing M particles 
and K particle-hole domain walls is realised in the steady state of the ASEP. Such 
configurations support K units of flux in the bulk, and so limj\r^oo Pn{Np, NJ) gives 
the joint probability of observing a density p and current J in the thermodynamic limit. 
To make use of earlier results for the ASEP, we note that 



xD yE 
1-xDl-yE 



^^^^"^ {D' + E'f. (3.95) 



:D'E' 



[1-x-y) 



2 



[1-x-y) 



2K " 



Then, from (3.93) and the relations (2.39)-(2.41) we have 

- ^(^"-y— (3.96) 
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in which is the normahsation for the ASEP as given in (3.38) but with the boundary 
rates a and (3 replaced by their primed versions 

a' = (3.97) 
1 — X — y 

P' = /""^ (3.98) 
1-x-y 

which arise from the requirement that D'\V) = {l/p')\V) and {W\E' = {l/a'){W\ in 
order for (3.96) to be correct. This expression is in agreement with Equation (4) in [104] 
which was obtained using almost the same approach. 

To analyse the joint current-density distribution Pn{M, K) in the thermodynamic 
limit, requires, once again, a careful analysis of the singularities in the generating 
function [104]. One anticipates expressions of the form P]\f{M, K) ^ exp[— A^r(p, J)] for 
large N, with the function r(p, J) = — limTv^oo ln-P/v(pA^, JN)/N taking its minimum 
value of zero at the thermodynamic values of the current and density given in Table 1. 
In the maximal-current phase one finds (up to additive constants relating to the 
normalisation of P) 

rip, J) = 2 Jin J + (p - J) ln(p - J) + (1 - p - J) ln(l - p - J) (3.99) 

and in the low-density phase 

r(p, J) = 2 Jin J + (p - J) ln(p - J) + (1 - p - J) ln(l - p - J) 

+ (l-p)lna + pln(l-a)-plnp-(l-p)ln(l-p) , (3.100) 

from which the result for the high-density phase follows from the particle-hole symmetry. 
In each case, the thermodynamic values give the desired minimum, and one notes the 
presence of non-Gaussian fluctuations. In [104], such fluctuations are suggested to be a 
signature of the long-range correlations that are present in the model. 



4. ASEP on a ring with two peirticle species 

In Section 2.4 we reviewed second-class particles and their various roles, and introduced 
a generalisation to two species of particles (see Figure 21): 

10^01 20A02 12^21. (4.1) 

The case a — (3 — 1 corresponds to species 1 being first-class particles and species 2 
being second-class particles. Let us stress here that we denote by a two-species exclusion 
process a system with two-species of particles and in addition vacancies. In some works it 
has been convenient to consider the vacancies as a species of particles and thus sometimes 
the two species system we study here has been referred to as a three species system in 
the literature. 



Matrix product expressions 

A microscopic conflguration is completely specifled by the set of occupation numbers 
where Tj = indicates that site i of the ring (relative to an arbitrary origin) is vacant 
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Figure 21. Asymmetric exclusion process on a ring with a second species of particle 
(shown as pentagons) . The labels indicate the rates at which the various particle moves 
can occur. 



and Tj = 1, 2 indicate occupancy by a particle of species 1 or species 2 respectively. We 
shall again deal with statistical weights (unnormalised probabilities) as in (2.34). On 
this occasion, 

F[Ti,T2, . . . ,TNj = (4.2) 

^N,M,P 

where here the normalisation is the sum of the weights over all configurations with M 
species 1 and P species 2 particles on the A^-site ring, since these quantities are all 
conserved by the dynamics. We extract a scalar from the matrix product appearing in 
(2.35) by the trace operation 

fin, r2, . . . , tn) = ti{X,,Xr, ■■■Xr^) (4.3) 

since this reflects the rotational invariance of the system. As before we use the notation 
D and E for the matrices relating to species 1 particles and vacant sites, and introduce 
A for species 2 particles. 

The reduction relations for these matrices are 

DE = D + E (4.4) 

AE = -A (4.5) 

a 

DA = ^A. (4.6) 

One notes the similarity with the relations (2.39)-(2.41) for the ASEP with open 
boundaries. Indeed, two sets of relations are equivalent if one chooses for A the outer 
product, or projector, 

A=\V){W\. (4.7) 

In particular, if one has a single species 2 particle [P = 1), the weights become 

/(n, T2,..., r^_i) = {W\X^,Xr, ■ ■ ■ Xr^^^lV) (4.8) 

if one chooses the location of the origin as being that of the second-class particle. 
Although the weight is identical to the matrix-product expression (2.35) for the case of 
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open boundaries, the distributions (2.34) and (4.2) are, in principle, different because 

the ensembles, and hence normalisations, are not the same. We shall return to this 
important point again below in Subsections 4.3 and 4.3.1. First, we show that this set 
of matrix product relations satisfies the stationary conditions for the system (4.1) (i.e., 
that illustrated in Figure 21). 

4-2. Proof of the reduction relations 

We use a simple generalisation of the algebraic approach of Sec. 3.1.2 to prove the 
reduction relations. The master equation can be written as a sum of terms, each relating 
to a pair of sites, as 

d ^ 

^/(ti, T2, . . . , Tat) = ^ hi^i+if(ri, T2, . . . , Tn) (4.9) 

1=1 

where in terms of matrix product expressions (4.3) for the weights 

hi^i+i tr(- • • [X^,X-r,+i]X^i+2 • • •) = • • ^n-i [W{Ti+iTi TiTi+i)Xi+iXi 

- w{nn+, ^ n+^n)x,,x,^jXr^^, ■■■), (4.io) 

in which W{t, t' — > r', r) gives the rate at which a neighbouring pair of sites in 
configuration (r, r') exchange. For two species of particles on a ring we have 

iy(10^01) = l (4.11) 
Vr(20 ^ 02) = a (4.12) 
1^(12^21) = /? (4.13) 

and all other rates zero. 

Introducing auxiliary matrices, Xt, r = 0, 1, 2, to generate a pairwise cancellation, 
we arrive at a sufficient condition 

w{n+in nTi+i)x^.^^x^. - winn+i Ti+iTi)x^,x^.^^ = x^.x^.^^ - x^.x^.^^ .(4.14) 

With a scalar choice for the auxiliaries and the set of rates (4.11)-(4.13) the 9 
equations in (4.14) reduce to three equations 

DE ^ED-DE (4.15) 
aAE ^EA-AE (4.16) 
I3DA ^AD-DA. (4.17) 

In this case, the choice D — —1, E —1 and A = yields (4.4)-(4.6) as required. 
4.3. Normalisation for a single defect particle 

On a periodic system the numbers of each species of particle are fixed, which puts 
constraints on the set of allowed configurations. For example, when P = 1, so that 
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there is a single species 2 particle (which we refer to as a defect), the normalisation for 
an A'"-site ring populated by M particles is given by 

11 1 N-l 

Zn,m=J2J2--- E {W\1[[t,E+{1-t,)DY^\V)Sj:^,^,m (4.18) 

ri=0T2=0 Tjv-i=0 j=l 

where Sij is the usual Kronecker 5-symbol and enforces the fixed particle number 
constraint. Although this expression can be evaluated directly, the general expression 
is rather cumbersome and the derivation relies on choosing a particular representation 
of the matrices [90, 105]. We consider instead the generating function 



(4.19) 



M=0 



which, by virtue of being now rather similar to (2.36) can be readily evaluated. To do 
this, one transforms to a new pair of matrices D' and E' via 



D^—{D'-1.) + 1 and E ^ y/^{E' - 1) + 1 



(4.20) 



which, by (4.4), obey the relations D'E' = D' + E', D'\V) = l/(3'\V), {W\E' = l/a'{W\ 
with the boundary parameters 



a' 



1 + ^(1-1 



and - 



lWH^-1 



(4.21) 



This transformation is very similar to that described in Sec. 3.7 to calculate the joint 
current-density distribution for the ASEP with open boundaries. Then, (4.19) may be 
written 

7V-1 



Fn{u) = {W\ (^^^ [D' + E'] + [^^ - 1] ' j \V) . 



(4.22) 



We can now use the formal power series approach of Section 3.2.3 to evaluate the 
generating function of Fn{u). Consider the formal series 



J^{u, z)^Yl ^'"''Fr.iu) = J2 {Wlz'^iV^C + iV^- 1)Y\V) 



N=l 



N=0 



;{W\ 



zu 



1/2 



1-^(01-1)' 



\v), 



where C = D' + E' . Now, we may write 



1/2 



1 - 



where 



Therefore, 



00(1 — cu) 



zwu 



i-z{^-iy ■ 



zu 



1/2 



1-^(01-1)2 



C 



= [1 - uj{u, z)E']''^ [1 - uj{u, z)D'] 
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(4.23) 



(4.24) 



(4.25) 



(4.26) 



and we deduce 

z) = 



-1 

(4.27) 



4.3.1. The grand canonical ensemble An explicit expression for Zjv,m could be 
extracted from Fj^{u) (4.19) via, for example, the residue theorem 

Zn,m = {u'^jFr.iu) = f ^^^iv(^^) . (4.28) 

Equivalently we can work in the grand canonical ensemble. To proceed we view (4.19) as 
a partition function for a system where the number of first-class particles can fluctuate. 
Then, u is a, fugacity that controls the density of occupied sites. We remark that it is, 
in fact, possible to generate this grand canonical ensemble dynamically by introducing 
certain specific dynamical rules which do not conserve the particle number as will be 
described in Section 7.3.2. 

In the grand canonical ensemble the mean particle number p is given by 

(M) 1 d 

p= lim lim — m— InFjv(w). (4.29) 

From our understanding of the grand canonical ensemble in equilibrium statistical 
physics, we anticipate that the both the mean number of particles (M) and its 
variance (M^) — (M)^ in the ensemble will be proportional to the system size N in the 
thermodynamic hmit. Thus, the particle number (M) will have fluctuations 0{N^^'^) 
about the mean and the density will be sharply peaked about p. The grand canonical 
and canonical ensembles should yield equivalent macroscopic properties. Since working 
at a fixed fugacity u is, in the thermodynamic limit, equivalent to working at fixed 
particle density p there is no need to invert the generating function (4.19). Instead 
(for large A^) one can substitute u for p via (4.29) in the expression for some physical 
quantity. 



4.. 4. Phase diagram for a single defect particle 

The phase diagram in the plane for the model with a single species 2 particle was 
determined in [90] by considering the exact expressions for the Zj^^m for N large. An 
alternative approach was presented in [84] where the Bethe ansatz was used to develop 
contour integral representations of Zn,m which were then evaluated in the large N limit. 
Here we will derive the phase digram from the asymptotic forms of the grand canonical 
partition function (4.19). We therefore proceed to evaluate the asymptotics of (4.19) 
from (4.27). 

We first note that the possible singularities of (4.27) as a function of z are (since u 
is positive) a square root at 

^0 = (1 + V^)"' , (4.30) 
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if a' < 1/2 a pole at u!{zo,) — o! and if /3' < 1/2 a pole at Loiz^ — (5' . The locations of 
the two poles reduce to 

ail — a) 



a{u - 1) + 1 



(4.31) 
(4.32) 



There is also an (irrelevant) pole at Zi = {1 — y^)~^. 

Using the results of Appendix B one can compute the following asymptotic 
contributions to Fn{u) : 

1 1 2a' 2(3' 



from 



yN 



from Zn 



from 



27r(l - Zo/ziY/^ a'-ip'-l 
1 [{1 — a)/a — ua/{l — a)] 
^ [(1 - a)/a - uil - (3)/(3\ 
1 [13 1 {I -P)- u{l - 
zg [(1 - a)/a - u{l - (3)/ (3] ' 

N 



2a' 



+ 



2(3' 



a' -I (3' -I 



(4.33) 
(4.34) 
(4.35) 



Since the asymptotic form of Fj^ is always ~ z , with z the dominant singularity, 



one can associate via (4.29) a density with each of the singularities 

d In zq v}-!'^ 



Pa = 



P/3 = 



u 



u 



u- 



du 
dlnza 

du 
d In 



1 + IiV2 

ua 

au + 1 — a 
u{l - (3) 



(4.36) 
(4.37) 
(4.38) 



du u{l -(3) +(3 

Note that since z^, < Zq and Zi3 < Zq a. contribution to F^^u) from Za or zp, if it exists, 
will dominate the contribution from zq. For example, the condition for the contribution 
from Za to exist is u^^"^ < {1 — a)/a. This becomes, using p — pa, (4.37) to give it as a 
function of p, 1 — a > p. Similarly, the condition for the contribution from zp to exist 
is w^/^ > (3/{l — (3) and this becomes using (4.38) (3 > p. 

Finally, we note that in both the z^ and zp contributions to Fn{u), Eqs. (4.34) and 
(4.35), there is a pole at 

u ^ /^(l-^) (4 39) 

{l-(3)a- ^ ' 

The presence of this pole means that as the parameters a and (3 change, u cannot cross 
Up. For the contribution coming from Za, one finds that as (3 decreases to p, u increases 
to Up. Similarly, for the contribution coming from Zj3, as a decreases to l — p, u decreases 
to Up. Therefore in the region a < 1 — p, P < p, u is fixed at m = Up. Actually, this holds 
strictly in the thermodynamic limit; to correctly achieve the density (p) = p one should 
choose u — Up + 7(0;, (3, p)/N, where the function 7 is some function to be determined, 
then let A?" — > 00. 

In the thermodynamic limit, we can can identify four phases as follows 
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Region 


Phase 


V 


P+ 


P- 


(i) 


PL 


l-2p 


l-{l-pf/a 




(ii) 


ER 


a — p 


P 


p(l - a) 113 


(iii) 


EL 


1-13-p 




P 


(iv) 


PC 


a — 13 


/3 


1 — a 



Table 2. Velocity, density in front of the defect and density behind the defect in the 
various phases of the single defect system 



(i) 1 — a < p < (3 z — zq — [1 — pY and u — — — ^ ^2 



(ii) 1 — a> p and p < 13 z = Za = Oi{l — p) and u 



Pf 

p{l - a) 



(1 - 

(iii) 1 — a < p and p > (3 z — z^ = (1 — I3){1 — p) and u 



[l-p){l-(3) 

(iv) 1 — a> p and p > /3 z — z^ — zp — a{l — f3) and u — — ^ 

We now calculate the velocity of the defect particle. In the canonical ensemble this 
is given by a(l — p+), where p+ is the probability of a hole in front of the defect minus 
(3p-, where p_ is the probability of a particle behind the defect particle: using (4.5) and 
(4.4) 



(4.40) 

(4.41) 
(4.42) 



Zn,m 

In the grand canonical ensemble the corresponding expression is 

a(a\E(uD + E)^-^\P) - P{a\(uD + E)^-^uD\(3) 

_ Fn-i{u) - uFn-i{u) 
Fn{u) 

We find using the form of the asymptotic behaviour, Fisr{u) ~ z~^ , that 

V z{l — u) for large . (4.43) 

The different phases are characterised by the values in Table 2. The phase diagram 
which results is shown in Figure 22, and the corresponding density profiles in Figure 23. 
Most of the phases are recognisable from the ASEP with open boundaries. First, the 
region (3 > p > 1 — a, which includes the case of a second-class particle {a — (3 — 1), 
has V — 1 — 2p and is marked PL (power law) in Figure 22. The effect of the defect is 
to cause a disturbance which results in an increase in the density in front of the defect 
and a decrease behind. The density profile as seen from the defect particle, decays as a 
power law with distance from the defect towards its asymptotic value p: 

(*-i^))"\ (4.44) 
^(i^y'^ . (4.45) 
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EL^ ' 




PC 


ER 




ER"^ 




1 . &- 



1-0 a 



Figure 22. Phase diagram for the ASEP with a single second-class particle on a ring. 

We note that this is essentially the same profile as that obtained for the ASEP with 
open boundaries in the maximal current phase. 

Since the bulk density far away from the defect is fixed at p in the periodic system, 
we do not have low- and high-density phases, as such, in this case. However, the 
characteristic density profiles of these phases, in which one has exponential decays to 
the left and right of the defect particle respectively, are present here in the phases marked 
EL^ and ER^ respectively. In the phase ER phase the defect only causes a disturbance 
in front of it and the density profile decays exponentially towards the asymptotic value. 
Similarly, in the phase EL phase the defect only causes a disturbance behind it. In the 
figure the + or — superscript indicates whether the density increases or decreases as one 
goes clockwise around the ring. The dashed line separating these regions with positive 
and negative density gradients has a + j3 = 1 and, in common with the ASEP with open 
boundaries, along this line the density profile is completely flat. 

The main novelty arising from the particle-number conservation is the PC (Phase 
Coexistence) phase, where a < 1 — p, /3 < p. Here, there is coexistence between a region 
of low-density (3 in front of the defect particle and high-density 1 — a behind the defect 
particle. At the point where the two regions meet there is a shock in the density. The 
shock is localised at a position a distance xN in front of the defect particle, where x is 
given by p = /3a; + (1 — a)x [90, 106]. The interesting point is that the shock exists in an 
entire region of the a-(3 plane as opposed to the open boundary condition case where 
phase coexistence only occurs along the line a = /? < 1/2. 

To understand the PC phase from the analysis, we note that there are two 
competing contributions to F^iu) from the pole in (4.27) at and the pole at zp. 
Using u = Up and (4.37,4.38) these correspond to densities (3 and 1 — a respectively. 
Finally we remark that in all phases the current of first-class particles J = p(l — p) 
except the PC phase where J = p{a — (3) + (3(1 — a). 

Finally, we mention that other defect particle dynamics have been studied [107-112] 
and we will discuss these models further in Section 7. 
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Figure 23. Density profiles as seen from the defect particle (shown as a vertical dotted 

line) in the four phases of the ASEP on a ring with a defect particle: power law (PL), 
phase coexistence (PC) , exponential decay to the right (ER) and left (EL) of the defect 
particle. 



4.5. Some properties of the second-class particle case 

The second-class particle case a = j3 = 1 was subject of a detailed study by Derrida, 
Janowsky, Lebowitz and Speer [78]. In addition to the case of a single second-class 
particle, results were obtained for several second-class particles and a finite density of 
second-class particles (as far as we are aware these computations have not yet been 
extended to the case of general a and (5). Here we briefly summarise some of their 
findings. 

For the case of two second-class particles the probabihty of finding the particles a 
distance r apart was computed and shown to decay as a 1 / r^/^ power law. Thus there 
is an effective attractive interaction and the two particles form a weak bound-state 
in which, despite the attractive interaction, the mean separation of the two particles 
diverges. The computations of [78] also considered a finite density, p2 of second-class 
particles and showed that in the limit p2 ^ ^ the environment around a second-class 
particle is different from case of a single second-class particle. Although in this hmit 
the global density is zero there must be a clustering of the second-class particles into a 
bound state. 

A finite system with first and second-class particles cannot support two regions of 
different densities. However via an astute construction [78,81] one can study shocks on 
such a system. The idea is that a system with finite densities pi, p2 of first and second- 
class particles can be viewed as a one-species system in two different ways. Treating 
the second-class particles as holes gives a TASEP with density pi whereas ignoring the 
difference between first and second-class particles, so that exchanges between them are 
immaterial, yields a TASEP with density pi-\- p2- Then one can construct shock profiles 
as follows: consider the density profile as seen from some second-class particle where 
in front of this second-class particle both first and second-class particles are treated 
as particles and behind the second-class particle the second-class particles are treated 
as holes. Then letting the system size go to infinity yields a profile where pi tends to 
pi -\- p2 for i large and p_j tends to pi for i large. Using this approach shock profiles were 
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Figure 24. The full parameter space of the partially asymmetric exclusion process 
with open boundaries. Particles can hop in cither direction on the lattice, and can 
enter and leave the system at both ends, all at different rates. 

obtained. Later these profiles were recovered from a calculation directly on the infinite 
system (see Section 7.2.3). 

5. The partially asymmetric exclusion process with open boundaries 

In this section we consider a generalisation of the ASEP to the situation where particles 
can hop both to the left and to the right, with different rates. We refer to this variant 
as the partially asymmetric exclusion process (PASEP). The dynamics are as shown in 
Figure 24, and one has particles on the lattice hopping to the left and right with rates 
p and q respectively, as long as the hard-core exclusion constraint is maintained at all 
times. Furthermore, particles enter the system at the left boundary at rate a, as before, 
but can also leave from the left boundary at rate 7. Likewise, particles leave the right 
boundary at rate /3, as before, but also enter there at rate 6. We shall once again define 
the unit of time by setting p = 1. 

Physical interest in this more general model comes from at least two quarters. First, 
there is a mapping between exclusion models and certain nonequilibrium surface growth 
models [113] (note, these are distinct from the equilibrium surface models discussed in 
Sec. 3.6). At large lengthscales, the shape of the interface is governed by a stochastic 
partial differential equation, the KPZ equation (see e.g. [31] for a review), that has a 
nonlinear term proportional to 1 — g. The presence or otherwise of this nonlinearity 
governs in turn the universal properties of the interface. Thus, by taking g — >■ 1 in the 
PASEP we should learn something of the crossover between universality classes [114]. 
The second source of interest lies in the special case where 7 = 5 = and g > 1. Then, 
the bias in the random walks performed by the particles opposes the current imposed 
by the boundary conditions. This gives rise to a new reverse-bias phase not seen in the 
ASEP. Furthermore, the matrix algebra that is used to solve the PASEP also arises in 
the context of ballistic reaction dynamics and has been used to calculate the late-time 
density decay in the system [115]. 

In Section 3.2 we discussed the three main approaches to calculation with matrix 
products. For the PASEP with open boundaries the most natural approach turns out 
to be diagonahsation of an explicit representation. Indeed, to our knowledge, results 
have thus far been obtained only through this approach [116-120]. 



59 



5.1. Quadratic algebra 

As for the ASEP, the stationary distribution is given through products of two matrices D 
and E that represent particles and holes respectively. Again, the product is contracted 
with a pair of vectors {W\ and \ V) to obtain the desired scalar weight, as in (2.35). The 
relations that allow the matrix products to be reduced are 

DE-qED ^D + E (5.1) 
{PD-5E)\V) ^\V) (5.2) 
{W\{aE--fD) = {W\ . (5.3) 

These can be proved to solve the stationary master equation either by considering terms 
that cancel between domains, as in Sec. 3.1, or using a slight modification of the purely 
algebraic approach of Sec. 3.1.2. We shall run over the latter briefiy here. 

As in Sec. 3.1.2, we decompose the terms in the master equation into a sum over 
terms relating to interactions at the boundaries, and between pairs of sites in the bulk. 
That is, 

7V-1 




where hi^i+i generate the terms arising from particle hops in the bulk, as in (3.9), and Kl 
and fiR perform the corresponding duties for interactions at the left and right boundaries 
as in (3.10,3.11). 

As in Sec. 4.2 we seek to telescope this sum, and we do so by assuming that the 
terms in the master equation generated by the h operators can be written using auxiliary 
matrices as 

ki+^{W\ ■ ■ ■Xr.X,,^, ...\V)^{W\--- [X,,Xr,^, - Xr,Xr,J ■.■]¥) (5.5) 
hL{W\Xr, ■■■\V)^ -{W\Xr, ■■■\V) (5.6) 
hn{W\..-X,jV)^{W\-..X,jV) . (5.7) 

Matching these right-hand sides up with the terms generated by the respective operators 
in the master equation yields 

= DD - DD = EE - EE (5.8) 

DE - qED =ED-ED = -DE + DE (5.9) 

{W\{-fD -aE)= - {W\E = {W\D (5.10) 

{SE-(3D) ^D\V)^-E\V). (5.11) 

Again, the scalar choices D = —1, E = 1 yield the consistent set of relations (5.1)-(5.3). 

It turns out that mathematical structure of the solution is much the same for the 
case with nonzero 7 and 6 as that with 7 = 5 = 0. However, in the former case the 
various boundary rates combine with one another in a complicated way, so we shall take 
7 = 5 = in the analysis that follows for ease of presentation. At the end of this section 
we shall explain how nonzero 7 and 5 were catered for in the works of [119, 120]. 
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5.2. Calculation of using an explicit representation 



As we observed in Section 2.3.1 of the Introduction, an argument was presented in [58] 
for the D and E matrices to have, in general, infinite dimension, at least for the case 
q — 0. This is also the case for general q, and one can find a number of different 
representations that satisfy the relations (5.1)-(5.3). The first of these that we shall 
consider has 



D 



1-q 



1 + 6 













1 + bq 













1 + 













1 + bq^ 



\ 



( 1 



E 



Co 1 + 










aq 
cT 









aq 







1 + ag^ 



V 



where the parameters a,h,Cn are given by 



a — 



Q 



a 



1, 6 = 



1-g 



and 



c„ = (l-g"+i)(l-a6g"). 
In this case, the appropriate choice for the boundary vectors is 

{Wi\ = (0| and l^i) = |0) 
where we have introduced a set of basis vectors 

|n) = (0,---,0,l,0,0,---,f 



(5.12) 



(5.13) 



(5.14) 

(5.15) 
(5.16) 

(5.17) 



for n = 0, 1, 2, . . . and their corresponding duals (n| by transposition. One can check 
exphcitly that the reduction relations (5.1)-(5.3) are satisfied. 

One way to evaluate the normalisation using this representation is to diagonalise 
the matrix C — D + E, as was done in Section 3.2.2 for the case q — 0. In analogy with 
(3.43) we parametrise the eigenvalues of C by 2(1 + cos^)/(l — q) where < ^ < 27r so 
that 



C|cos^) = 



2(1 + cos ^) 



cos 9) . 



Writing 



cos 6') = ^/„(cos6l)|n) 

n=0 



(5.18) 



(5.19) 
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implies from (5.12) and (5.13) that 
2(1 + cos e) 



1 



Putting 



-fn{cOs9) 



(5.20) 



Y^{v^/n-i(cose) + [2 + (a + b)q'']U{cose) + v^^/n+i (cos ^) } . 



f„{cos9) 



PnjcOSe) 

nn—l I — 
A;=0 V '^fc 



(5.21) 



this three-term recursion simphfies to 

2cos^P„(cos^) = (l-g")(l-a6g"-i)P„_i(cos^) + (a+6)g"P„(cos^)+P„+i(cos^) .(5.22) 

Given the boundary conditions P_i(cos^^) = 0,Po(cos^^) = 1, this recursion defines the 
set of Al-Salam-Chihara polynomials [121]. The relevance of these polynomials to the 
PASEP was noticed by Sasamoto [116]. 

Although these polynomials are rather complicated, we need to know very little 
about them to evaluate the normahsation Zjsi — (l^|C^|y), other than that they satisfy 
an orthogonality relation [116, 121] 



1 {(i,uh\(i) 



2i0 „-2/(?. 



d^ 



27r (g, ab; q)n Jo (oe^^ ae'^^, be'^, be''^; q). 



-Pn{cOs9)Pm{cOs9) = Sn,m, 



(5.23) 



which holds when \a\ < l,\b\ < 1 and \q\ < 1. Here we have invoked the notation 



n-l 



(oi, 02, ... , a„; q)n = JJ(1 - aiq''){l - a2q^) • • • (1 - a„g* 



(5.24) 



fe=0 



since the quantity on the right-hand side is ubiquitous in the theory of ^-series, of 
which the Al-Salam-Chihara polynomials are an example [121]. Using this orthogonahty 
relation, we can construct a representation of the identity 



1 



{Q,ab;q)c 
27r 



d^- 



(e2^^e-2^^;?)c 



cos ^) (cos 6* I 



(5.25) 



/o (ae^^,ae-^^,&e^^,6e-'^;g)c 
where (cos^| is simply the transpose of | cos^). 

An integral representation for the normalisation Z is then obtained as follows. First, 
one incorporates this identity into the definition of Z via Zjv — {W\C'^'m.\V) . This yields 



Zn — 



{Q,ab]q)c 
2n 



d9- 



(iy|c^| COS ^) (cosily) . 



(5.26) 



We then use the fact that, by construction, C\ cos 6*) = 2(1 + cos^)/(l — q) to obtain 

N 



z ^ 



(g,Q^;g)c 

27r 



d^- 



ae'^, ae-'^ 6e*^ be''^; q)c 



.{W\cos9) ^ ^ (cos^lV^). 

L 1-9 J 

(5.27) 

Finally we note that (W^|cos^) = (O|cos^) = 1 and likewise that (cos^|\^) — 1. We 
find then an integral representation of the normalisation 
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2(1 + cos( 
1-9 



N 



.(5.28) 



In the limit g — > this can be written as 

= izi^ r dd- ^- [2(1 + cos e)f (5.29) 

TT (l + a2-2acos^)(l + 62-26cos^) ^ ^ ^ ' 

which can be shown to agree with the finite sum (3.38) as it should — see Appendix C. 
In [117] it was also shown that (5.28) is equivalent to a finite sum expression for general 
q. However, the resulting expression is rather cumbersome, and so in Section 5.3 we will 
concentrate on the procedure to extract the asymptotics of the normalisation directly 
from the integral representation (5.28). 

5.2.1. Relation to Motzkin Paths The representation given as (5.12) and (5.13) has 
an interpretation in terms of random walks and surfaces, similar to that described in 
Section 3.6. To see this, let the vector \n) represent a height n above some origin. A 
matrix element (m|C/|n), where C/ is a product oi £ D and E matrices can then be 
interpreted as the weight of a set of paths of horizontal length i connecting a point at 
height m (at the left-hand end of the surface) to a point at height n (at the right-hand 
end of the surface). The particular choice of boundary vectors {Wi\ = (0| and \Vi) = |0) 
thus corresponds to a set of paths that begin and end at the origin. 

The bidiagonal forms of the D and E matrices imply that a step of the surface 
must be fiat, a diagonal up-step or a diagonal down-step. Furthermore, the fact that 
the matrices are semi-infinite means that the paths may never go below the origin. 
Finally, the matrix product Zn — (0|(-D -|- £')-^|0) describes the ensemble of all possible 
such paths of length N, and has up- and down-step pairs connecting height n and n + 1 
weighted by c„, and two different types (or "colours") of horizontal segments, coming in 
with weights 1 + ag" and 1 + bq^, where again n is the height of the segment above the 
origin: see Figure 25. These paths sometimes appear in the combinatorial literature as 
bicoloured Motzkin paths. 

We remark that each bicoloured Motzkin path can be mapped uniquely to a Dyck 
path, which has only diagonal up- and down-steps as described in Section 3.6. This is 
achieved by mapping each segment of a Motzkin path to two consecutive segments of 
a Dyck path. Specifically, diagonal up- and down-steps are mapped to a pair of up- 
and down-steps respectively, whilst the two colours of horizontal steps are mapped to 
up-down and down-up pairs. In order that the resulting path contacts the origin only 
at its two ends (i.e., is a Dyck path) two additional up-steps are needed at the left end, 
and similarly two down-steps at the right-sec Figure 25. Thus a Motzkin path of length 
N corresponds to a Dyck path of length 2(A^ + 2). Although a similar treatment of 
the PASEP thermodynamics to that described for the ASEP in Section 3.6 should in 
principle be possible using this approach, to our knowledge it has not yet been achieved. 
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(a) 



(b) 



■-Or 



\ 



/\/\ 



Figure 25. (a) A bicolourcd Motzkin path. Each step may increase or decrease 
the height of the surface by one unit, or leave it unchanged. These latter horizontal 
segments come in two colours, denoted here by the flat and crenellated segments, (b) 
The same path mapped onto a Dyck path, wherein each Motzkin segment becomes 
two Dyck segments. Two additional segments (shown dotted) are needed at both ends 
to prevent the Dyck path from crossing the origin. 



5.3. Asymptotics of the normalisation from the integral representation 

To analyse the asymptotics of the normahsation from its integral representation (5.28) 
we first rewrite it as the contour integral 



{q,ab;q)oo f dz {z'^,z ^;g)oo ^ 2 + z + z ^"i^ 

1-q 



Am z {az,a/ z,bz,b/ z]q)oo \_ 1 — ? J 

where K is the positively oriented unit circle centred on the origin. Recall that the 
orthogonality relation (5.23) was valid only when \a\ < 1 and \h\ < 1. In this case 
the contour K encloses two sequences of poles at z = a, ga, q%, . . . and z = b, qb, q^, ■ ■ ■ 
along with a number of singularities at the origin, as shown in Figure 26 (i). We can 
continue the function Zn to a > 1 and/or 6 > 1 by ensuring the contour K is deformed 
still to enclose these poles, and further the exclude those at z — l/a,q/a, q^/a, . . . and 
z — 1/h, q/b, q^ /b, ... as shown in Figure 26(ii). 

In the former case, where \a\ < 1 and |6| < 1, the contour K can be continuously 
deformed without passing through any poles to coincide with the line that passes parallel 
to the imaginary axis through the saddle point in the integrand at 2; = 1. See Figure26(i). 
Then the normalisation is well approximated for large N by the expression obtained from 
applying the saddle-point method. This is [116,117] 

4 {q;q)l{ab;qU 4^ ..^.s 

Let us now consider the case a>l>ga, |6|<1 where the deformed contour 
shown in Figure 26 (ii) cannot be continuously deformed into that used in the saddle- 
point method. One must therefore correct the result by applying the residue theorem to 
the poles at 2; = a and 2; = 1/a so that the former is included in the integration contour 
(as required) and the latter excluded. This procedure gives a correction 

ZN-^-llP^i2 + a + l/a)'' (5.32) 

which, for sufficiently large A^, overwhelms the contribution from saddle-point expression 
(5.31) because 2 + a+ l/a > 4 when a > 1. Therefore, the normalisation is well- 
approximated by the contribution (5.32) from the pair of poles at 2; = a and z = 1/a. 
As a is further increased, more and more poles need — in principle — to be accounted for 
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Figure 26. Integration contour K for the representation (5.30). Tlie x symbols show 
the first of each of the sequences of poles 2: = a, ga, g^a, . . . , and z — a. qb, q'', . . . that 
are required to be inside K, and the + symbols those that must be outside, (i) When 
\a\ < 1 and |5| < 1, K can be continuously deformed into the contour for the saddle- 
point integration (the vertical line) as shown, (ii) When \a\ > 1, this deformation is 
no longer possible, and it is necessary to add to the saddle-point result the difference 
between the residues at z = a and z — I /a. 



Region 


Phase 


Current J 


a < l3,a < ^ 
/? < < 


Low-density (LD) 
High-density (HD) 
Maximal current (MC) 


a{l—q—a) 
1-g 

l-q 
l-q 
4 



Table 3. Asymptotic currents in the various phases of the PASEP. 



in the deformation of the contour. However, one can verify that all further corrections 
grow less fast than {2 + a+ and so (5.32) remains the leading contribution to the 

normalisation. 

When \a\ < 1 but 6 > 1, we arrive at the same expression (5.32) for the 
normalisation but with a and b exchanged. When both a and b are larger than 1, 
both of these contributions are present. Here, (5.32) dominates when a > b, and the 
corresponding expression with a ^ b when a < b. Thus, for the PASEP we have the 
same three phases as for the ASEP. In each of these the current can be determined using 
the matrix expression 

J ^ {W\C'-\DE~qED)C''-^-'\V) ^ Zn-1 
Zn Zn 
where the reduction relation (5.1) has been used to yield the same ratio of normalisations 
that gave the current for the ASEP. Rewriting the inequalities on a and b given above in 
terms of a and f3 we find the phase diagram shown in Figure 27 with the corresponding 
currents given in Table 3. Notice the similarities with Figure 5 and Table 1. 

5.4- Density profiles from the q-deformed harmonic oscillator algebra 

Whilst the matrix representation given in (5.12) and (5.13) provides a quick route to 
the normalisation (5.28), it is less suited to the task of calculating density profiles. 
Here, progress is made by using a representation that involves the g-deformed harmonic 
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l-q 
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Figure 27. Phases of the PASEP over which the particle current is analytic. 



oscillator algebra, a well-studied mathematical object for which many results are known 
[122,123]. 

Central to this approach is a pair of creation and annihilation operators and d 
that satisfy a g-commutation relation 



dd^ — gd^d = 1 — q 
and act on basis vectors according to 

dV) = a/1 - + 1) 



a\n) = y^l — q^\n — 1) . 



If we take 



D 
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'1 + d) and E 



1 



(l + dt) 



(5.34) 

(5.35) 
(5.36) 

(5.37) 



1 — g ^ — Q 

we find from (5.34) that the reduction relation (5.1) is satisfied. Thus this representation 
looks like 

/ 1 

1 v^l -g2 

1 v^l -g3 

1 



D 



1 



\ 



(5.38) 



E 



l-q 



/ 1 ■■■ \ 

v^r^ 1 

v^l-g2 1 

v^l-g3 1 



(5.39) 



In order that (5.2) and (5.3) are also respected (bearing in mind 7 = 5 = 0) we require 



{W\n) = K 



—j^=^= and {n\V) = k = 



(5.40) 
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where the parameters a and h are as previously given in (5.14). Thus this representation 
is the gcnerahsation to g of that used in Section 3.2.2. The constant k appearing 
in (5.40) is fixed by the convention that = 1- That is, we require 



00 



4 = E 7^ • (5-41) 

This sum, which converges when \ab\ < 1, features prominently in the literature on 
g-series as it is the g-analogue of the exponential function [121]. When it converges, it 
can be expressed as an infinite product and = {ah; q)oo- 

Diagonalisation of the matrix C = D + E proceeds in much the same way as 
described above in Subsection 5.2. Here, the eigenfunctions are a vector of functions 

|cose) = ^^^^=i|n) (5.42) 

with eigenvalue A(cos^^) = 2(1 + cos^^)/(l — g) as before. The polynomials Hn{cos9) are 
found to satisfy the three-term recurrence 

2cos9H^{cos9) = (1 - g")//„_i(cos ^) + Hn+i{cos9) n>0 (5.43) 

using the same approach as previously. Taking Hq = 1, allows the functions Hn{cos9) 
to be identified as g-Hcrmite polynomials, the properties of which are discussed in 
[7, 93, 121]. Of particular importance is the g-exponential generating function of 
these polynomials, since this allows computation of the scalar products {W\ cos 9) and 
(cos^|y) that appear during the calculation of the normahsation [116,117]. One finds 
for the former 

/TT.I /IV a"if„ (cos ^) 1 /r..N 

and the corresponding expression with a — > 6 for the latter. With the knowledge of 
the orthogonahty relation satisfied by the polynomials [93, 121] one can write down the 
identity 

j^^Moo rd^(e2ie^e-'*';g)oo|cos^^)(cos^^| (5.45) 
27r Jo 

which is the final part of the jigsaw that is an expression for the normalisation. One can 
verify that the same integral (5.28) is obtained using this alternative representation. 

The key benefit offered by this approach is that expressions for the difference in 
density between neighbouring pairs of sites 

Ai = pi- pi+i = ^ ^ ^ ^ (5.46) 

involve the matrix DC — CD — DE — ED which is diagonal in this representation. This 
one sees from (5.37) and the definitions (5.35) and (5.36). Specifically, 

, ,r^„ ^^1, > (n\\aa^ — a\\m) 

{n\ [DE - ED] \m) = ^ ^ ^ (5.47) 



1-g 

[(1 - - (1 - g")] 

1-g 
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Sn,m = q''Sn,m ■ (5.48) 




1+q 2 



Figure 28. Full phase diagram of the PASEP showing additional sub-phases that 
have different density decay forms. 



To express the density gradient (5.46) in terms of integrals over the g-Hermite 
polynomials, one inserts two of the identities (5.45), one in front of the combination 
[DC — ED] and one behind. This gives rise to a double integral which is analysed for 
large N and i in much the same way as discussed above in Subsection 5.3 but with 
additional complexity stemming from the fact that one has the deformation of two 
contours to contend with. 

The detailed analysis of this double integral was conducted in [118]. As with the 
ASEP, one finds additional sub-phases relating to different boundary profiles within the 
high-density and low-density phases. The full phase diagram is as shown in Figure 28 
and the density profiles one obtains summarised as follow. 

• Maximal-current phase a > (3 > Once again, one has the universal 
power-law density profile seen both in the ASEP with open boundaries and on a 
ring for sufficiently large a and f3. That is, near the right boundary one has the 
decay 

(5.49) 



2 V 0U 

and a similar decay from the left obtained by invoking the particle hole symmetry 
that was also present in the ASEP, viz, pi{a,P) = 1 — pN+i-i{P,Oi). 

• Low-density phases a < < /5. In all of the low-density phases there is a 

region of constant density p = < | that extends from the left boundary into the 
bulk. Near the right boundary the decay takes three different forms which define 
three sub-phases. 

— LD-I a < P < < Here there is an exponential decay with a 



characteristic length given by 

-/3{l-q-/3) 



C = In 



a{l — q — a) 



(5.50) 
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This phase is present in the totally symmetric case g = (see Figure 17). 
LD-II ^Yf^ < a < > This phase, which is also present when 

g = 0, has an exponential decay characterised by 



In 



4a{l — q — a) 



(5.51) 



which is additionally modulated by a power- law decay with exponent |, as in 
the totally asymmetric case. 



1+9 ■ 

purely exponential decay with a characteristic length 



LD-III /3 > ^,Q! < ■ When g > 0, a new phase appears that has a 



= In 



{q + a] 



(5.52) 



• High-density phases f3 < ^-^,/3 < a. The properties of the high-density phases 
are obtained from the particle-hole symmetry. These all have an exponential decay 
at the left boundary to a constant bulk density that extends to the right boundary. 

5.5. Reverse bias phase 

The foregoing analysis applies to the case g < 1, since this is a restriction on the 
representations of the identity (5.25) and (5.45) that were required in the analysis. As 
we previously noted, the case where g > 1 is of interest since then the particles in the 
bulk prefer to hop to the left, even though the boundary conditions enforce a current 
that flows in the opposite direction. 

Although the g-Hermite polynomials have an extension to g > 1 [124] it is unclear 
how to handle the resulting integral representation of the normalisation and other matrix 
product expressions [117]. An alternative approach, followed in [117], is to note that, if 
one were actually to evaluate the integral (5.28) by performing, say, a Laurent expansion, 
one would end up with a sum containing a finite number of terms. This is because this 
integral is another way to express the result of a direct reordering using the reduction 
relations (5.1)-(5.3), a procedure that is known to terminate. Such an expression is 
then valid for any g, not just g < 1 as was required to derive the integral representation 
in the first place. In other words, we can access the asymptotics for g > 1 by first fixing 
N at some finite value with g < 1, put g > 1 in the result and only then at the very end 
take N ^ oo. 

It is possible, using a clutch of identities from the theory of g-series, to find the finite 
sum implied by (5.28) [117]. Wc do not quote this formula — a complicated mixture of 
ordinary and g-combinatorial factors — here. Rather, we simply state that the leading 
term in the normalisation for g > 1 is, for large N [117], 



TV 



y A( m( 1 (g-l + a)(g-l + /?) \ 1^2 
Z^^A(g,a,/3)(^^y J g. (5.53) 

where A{q, a, (3) is a function that is independent of N. Note that this result holds for 
all values of a and (3: there is only one reverse-bias phase. For (5.33) we see that in this 
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phase, the current vanishes asymptotically as 

(5.54) 

\(q-l + a)(q-l + P)J ^ ' 

Although the density profile in the reverse bias phase has not been calculated 
explicitly, one can conjecture its form using a simple argument. First we expect most 
of the particles to be in a domain by the left boundary, and most of the vacant sites by 
the right. Now, if the interface between these particle- and hole-rich regions were closer 
to the right boundary than the left, a particle at the interface will escape from the right 
boundary in a shorter time than a hole from the left. Therefore, in this instance the 
interface will tend to drift to the left. On the other hand, if the interface were closer 
to the left boundary it would, by the same argument, drift to the right. Therefore, 
the lattice will typically be half full in the steady state, a feature noticed in computer 
simulations for small systems. 

5.6. Entry and exit at both boundaries 

We finally briefly mention the case of the PASEP where particles may enter and leave 
the system at both boundaries, as shown in Figure 24. A representation similar to 
that used in Subsection 5.2 was established in [119], where it was also shown that the 
functions appearing in the eigenvector of the matrix C are a complete set of Askey- 
Wilson polynomials [121, 125]. These are a set of very general polynomials that contain 
both the Al-Salam-Chihara and g-Hcrmitc polynomials encountered above as special 
cases. They have four parameters related to the transition rates a, /5, 7 and S as follows 

a± = 7^ [(1 - g - a 7) ± - q - a + -f)"^ + Aaj] (5.55) 
Za I J 



[1 - q - 13 + 6) ± - q - 13 + 6)^ + Apd . (5.56) 



As in Subsection 5.2, the boundary vectors are {W\ — (0| and \V) — |0), and the 
polynomials satisfy an orthogonality relation similar to (5.23), but with additional 
factors appearing to accommodate the new parameters [121, 125]. These carry through 
to the expression one eventually obtains for the normalisation 
^ _ {q, a+a-, h+h^, a+b+, a+b^,a-b+, a-b-] g)oo ^ 
^ 27r(a+6+a_6_, g)oo 
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a+e^^a+c-»'',6+e*^,6+e-*'',a_e^^a_e-*'',6_-e*^,6_e-'^;g)oo L ^ ~ Q J 
One recognises a structure similar to (5.28) and the asymptotic analysis via a contour 
integral detailed in [119] is basically the same as described in Subsection 5.3. When 
q < 1, one finds the three familiar phases: a maximal-current phase when a+ > 1, 6_|_ > 1; 
a low-density- phase when a+ > 1 and a+ > 6+; and a high-density phase when 6+ > 1 
and 6+ > a+. When g > 1 the same three phases are manifested as the model is invariant 
under the exchanges {a, (3, 7, 5, q) {q~^S, q"^^, q~^P, q~^Oi, q~^), at least when all four 
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boundary rates are nonzero. Since the physical properties of these phases are very 
similar to their counterparts in the ASEP with open boundaries, we do not discuss 
them further here. Further details, including an analysis of the bulk densities through 
the introduction of a fugacity as described in Subsection (4.3.1), are given in [119]. As 
further shown in [120] , the eigenvalue spectrum of the matrix C can be used to determine 
density profiles and correlation functions. These show once again the universal power- 
law decay in the maximal current phase, and various exponential decays in the high- 
and low-density phases. 

6. Macroscopic density profiles for open boundeiry ASEP 

In this review we have so far mostly concentrated on how the current and density profiles 
in driven diffusive systems are calculated from the matrix product expressions for the 
steady state. Since one has to hand the entire stationary distribution of microstates, it 
is in principle possible to calculate any macroscopic quantity of interest. For example, 
the existence of a matrix product solution has recently admitted the calculation of the 
Gibbs-Shannon entropy in the symmetric exclusion process with open boundaries [126]. 
In this section we take a closer look at a quantity that has received a great deal of recent 
attention in the context of open driven-diffusive systems, namely free energy functionals 
which characterise the probability that a deviation from the most likely density profile 
is seen in a macroscopically large (but finite) system. 

To be more precise, we consider macrostates corresponding to a density profile 
p{x) where x lies in the (one-dimensional) interval [0, 1] which contains some large 
number N of lattice sites. As we shall see explicitly below, each macroscopic profile 
p{x) can be realised in number of ways that grows exponentially with the number of 
microscopic degrees of freedom N . This exponential growth leads to the key property of 
statistical mechanical systems that in the thermodynamic limit, — oo, one particular 
macrostate — let us call this p{x) — is very much more likely to be realised than any 
other (except possibly at a phase transition). That is, the probability P[p] of seeing a 
particular macrostate behaves for large N as 

P[p] ~ exp (-Arjr[p]) (6.1) 

where T[p] is a functional of the density profile that vanishes when p{x) — p{x). In 
an equilibrium system, jF[p] is essentially the free energy of the macrostate, relative to 
the most likely profile. Of course, in a nonequilibrium steady state, the probability of 
seeing a deviation from the most like profile can be measured, and so by analogy with 
the equilibrium case, we refer to the object T[p\ that results as a free energy functional. 
Progress in calculating such free energy functionals and related quantities has been 
summarised in the lecture notes of Derrida [127]. 
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6.1. Free energy functional for the symmetric simple exclusion process (SSEP) 

In this subsection, we will sketch the derivation of the free energy functional for the 
nonequilibrium symmetric simple exclusion process (SSEP) with open boundaries that 
has been achieved through the matrix product solution for the steady-state distribution 
[128, 129]. This version of the model has symmetric hopping of a single species at unit 
rate in the bulk, entry at the left and right boundaries at rate a and 5 respectively, 
and exit at the left and right boundaries at rate j3 and 7. In other words, we have the 
dynamics illustrated in Figure 24 of Section 5 with the bias q = 1- 

In order to illustrate the procedure by which the free energy functional was first 
calculated in [128, 129] from the matrix product expressions, we shall follow it here for 
the much simpler case in which the model parameters are chosen so as to give rise to 
an equilibrium steady state. To find this set of parameters, we employ a Kolmogorov 
criterion [9, 130]. This states that if, for every loop Ci —> C2 Cn Ci in 

configuration space, the equality 

W{C^ ^ C2)W{C2 ^ C3) • • • W{Cn-l ^ Cn)W{Cn ^ Ci) = 

W{Cn ^ Cn-l)W{Cn-l ^ Cn-2) ' ' ' W {C2 ^ C^)W {C^ ^ C„) (6.2) 

holds, the steady state of the system is an equilibrium state for which the detailed 
balance relation 

f{C) _ WjC ^ C) 

f{C') W{C ^ C) ^ ' ' 

also holds for every pair of configurations C and C between which transitions occur at 
nonzero rate. Here we find that (6.2) is automatically satisfied unless, in total, some 
nonzero number m particles exits at one boundary and re-enters at the other to return 
to the starting configuration going one way round the loop. Then, the ratio of the two 
products in (6.2) is {a(3/jd)"^. Therefore, to realise an equilibrium steady state one 
must have = jS. 

Since the bulk hopping occurs at the same rate in both directions, the detailed 
balance relation (6.3) implies that all configurations that have M particles on the lattice 
are equally likely in the steady state. By considering a pair of configurations that differ 
by the addition or removal of a single particle at the left boundary, one finds that if /m 
is the steady-state weight of any one of the M-particle configurations, 

/m+1 = -/m . (6.4) 

7 

Hence, the probability of seeing a particular configuration ti,T2, . . . ,rN (where, as 
previously, = 1 if site i is occupied, and zero otherwise) is given by the product 

N /aNTi 

P{n,T2,...,T^)^llj^. (6.5) 

i=l 7 

Once the underlying distribution of microstates has been established, the free energy 
functional is obtained by taking an appropriate combination of thermodynamic and 
continuum limits. To this end, we divide the system into a number n of boxes, box j 
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containing Nj — i/jN sites and Mj — PjNj particles. The j pairs {yj,pj) can thus be 
used to specify a macrostate of the system in the hmit N ^ oo. To find the probabihty 
of this macrostate we note that 



P{N,,M^;N2,M2r 



f[Pbox(iVi,M,) 



(6.6) 



where Phox{N, M) is the probabihty that N particles are in a box of size M in the steady 
state. This is given by 

Pbox(iV,M)= f (6.7) 
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where the second expression, valid for large A^, has been obtained using Stirling's 
approximation. Inserting this into (6.6), and taking the thermodynamic limit, one 
finds that 
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(6.9) 



in terms of the intensive box sizes yj and densities pj. The continuum limit is now 
straightforward: we take yi = 5xi and j ^ oo where Sxi is some portion of the 
interval x G [0, 1] such that dxi = 1 and the density profile p{x) becomes a function 
of x. This procedure yields 

p{x 



dx 



p{x)ln^ + {l- p{x))\n- 



P 1 - P 

for the free energy functional, where we have introduced the parameter p = 
is easily verified that J-'[p{x)] vanishes if p{x) = p Vx e [0, 1], indicating that the flat 
profile is the most likely profile in the thermodynamic limit as one would expect at 
equilibrium. We also note that the free energy functional is local in the density field 
p{x); that is, one can express (6.10) as an integral over a free energy density f{x) that 
depends only on the local properties of the density field p{x) at the point x (and nowhere 
else) . 

When the equilibrium constraint on the boundary rates af3 = ^5 is relaxed, this 
local property of the free energy density no longer holds. The reason for this is that when 
the system is divided into boxes, the joint distribution of box sizes and occupancy does 
not factorise as in (6.6). That is, there are long-range correlations in the density profiles 
that are implied by the matrix-product structure of the steady-state weights. This 
makes evaluating the free energy functional by taking the appropriate combination of 
thermodynamic and continuum limits a much more complicated exercise. We therefore 
omit the technical details of the calculation here, referring the reader to [129] for a full 
account. 

The way in which the free energy functional is obtained in [128, 129] is to introduce 
a pair of fugacities for each box, one of which controls its mean size yi and the other 



(6.10) 
It 
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its density pi. This procedure is similar to the use of the grand canonical ensemble 
described in Section 4.3.1, albeit with a much larger number of fugacities, and yields 
the generating function of the joint distribution that extends to the nonequilibrium case 
the expression (6.6). The values of the fugacities are then fixed by finding the saddle- 
point of this generating function. It turns out that the resulting expression for the free 
energy functional, 



T[p] = dx \p{x) In 4^ + (1 - p{x)) In + In ^'^''^ 



a{x) cT(l)-a(0) 

involves a companion function a{x) which is related to the fugacities described above, 
and that is fixed by the saddle-point analysis as the function that maximises the 
functional !F[p] for a given profile p{x). One can verify, by asking for a vanishing 
functional derivative that this function must satisfy the nonlinear differential 

equation 



,(6.11) 



(j"{x] 



a{x) + a{x){l - = P{x) ■ (6.12) 

Additionally, the detailed analysis presented in [128, 129] shows that (t{x) must further 

be a monotonic function, and satisfy the boundary conditions 

a 5 
(7(0) = and ail) = . (6.13) 

We note that the monotonicity requirement on a{x) ensures that the argument of 
the final logarithm appearing in (6.11) is positive, and further that the value of the 
companion function (7{x) at some point x will in general depend on the entire profile 
p{x) through the differential equation (6.12) — the free energy density is non-local, as 
previously claimed. 

6.2. Effective local thermal equilibrium and additivity principle 

The companion function a{x) that appears in the nonequilibrium free energy functional 
(6.11) is somewhat mysterious at a first encounter. However, a physical meaning can 
be ascribed to this function, starting from the observation that the values imposed at 
the boundaries coincide with the mean densities of the left and right boundary sites in 
the steady state, and hence also the densities of the particle reservoirs at each end of 
the system. In the bulk, one finds that an infinitesimal segment of the system located 
at a point x G (0, 1) can be thought of as being in local equilibrium with a reservoir at 
density a{x). 

To see this, it is convenient first to generalise the functional (6.11) to a density 
profile p{x) defined on an arbitrary interval x e [a, b] in such a way that if the profile 
is translated and scaled so that it covers a different interval x e [a', b'], the free energy 
per unit length is unchanged. The functional that has this property and coincides with 
(6.11) for the case a = 0, 6 = 1 is 

.(6.14) 



J a 



p{x) In + (1 - p{x)) In + In ^ ' ^ 



a{x) 1 — a{x) a{b) — a{a) 
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In the limit of an infinitesimal interval, 6 — > a, this expression approaches the 
corresponding expression for the equilibrium system 

/}\ 'I- I I — /t\ I 

(6.15) 



J a 



p(x)ln^-V^ + (1 -p(x))ln ^ 



P 1 -p 

if we identify (T{a) with p, the density of the reservoirs coupled to the equilibrium system. 

The observation that an effective local thermal equilibrium applies does not itself 
allow one to recover the full expression (6.14) for the free energy functional since a 
priori one does not know what densities the effective intermediate reservoirs should 
have. However, it was noticed in [129] that one can construct the free energy functional 
through an additivity principle that involves the modified free energy 

^[a,6]([p];o"(a),o-(6)) = J^[a^b]{[p]; Pa, pb) + (6 - fl) lu J(cr(a) , a(6) ) (6.16) 

where J {a (a), (j{b)) is the steady-state current through a system of length b — a coupled 
to a reservoir of density a{a) at the left boundary, and of density a{b) at the right. For 
the SSEP, this current is 

J = ^^Z^ , (6.17) 
b — a 

where we have assumed that a{b) > a{a). (Since the dynamics are symmetric, the case 
a{a) > a{b) can be treated by making the replacement x — > 1 — x.) The additivity 
principle given in [129] relates the free energy functional for a system to that of two 
subsystems created by inserting a reservoir at a point y e (a, b). It reads 

^M([p];(^(a),(7(^)) = max{7^[„,j^]([p];a(a),cT(|/)) +H[2,,b]([p];cT(|/),cT(6))} . (6.18) 

That is, (j{y) is the density one should choose for the reservoir placed at the point y so 
that the combined free energy of the two subsystems is maximised. One can verify by 
recursively subdividing the subsystems created by inserting intermediate reservoirs, and 
by assuming that each infinitesimal segment created in this way is in a local equilibrium 
with its boundary reservoirs, that the expression (6.11) results. The fact that one is 
looking for a maximum in (6.18) further implies that the reservoir densities a{x) must 
satisfy the differential equation (6.12). Thus the additivity principle (6.18) in tandem 
with the effective local equilibrium property and the free energy functional given by 
(6.12) are equivalent. 

6.3. Properties and applications of the free energy functional for the SSEP 

It is worth reiterating a few properties exhibited by the free energy functional (6.11) 
noted in [129]. First, as is required, the equilibrium version (6.10) is recovered in the 
limit where the right boundary reservoir density a(l) approaches that of the left (t(0). 
To see this, we write (without loss of generality) for the effective reservoir density in the 
bulk 

a{x) = (j(0) + [(j(l) - a{Q)]x + 5a{x) . (6.19) 
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If one substitutes this expression into the differential equation (6.12) that governs (t{x), 
one finds that the function da{x) must be proportional to [o"(l) — o"(0)]^ if a{x) is to 
remain bounded as a{l) o'(O). In turn, this implies that a{x) a{0) = a(l) in 
this limit, yielding the desired equilibrium free energy functional, with p equal to both 
boundary reservoir densities. 

It is also straightforward to show that the most probable profile coincides with 
the stationary profile, which for the SSEP is simply the hnear function p{x) — 
(t(0) + [(7(1) — a{0)]x. This is achieved by asking for the functional derivative of (6.11) 
with respect to p{x) to vanish. That is. 



and hence we must have a{x) = p{x). The differential equation (6.12) then implies that 
the second derivative of a{x) must vanish, which (along with the boundary conditions) 
results in the optimal a{x) and p{x) coinciding with the linear stationary profile. It is 
easy to see that then (6.11) is zero, indicating that this profile appears with probability 
1 in the thermodynamic limit; it can also be shown [129] that any other macroscopic 
profile appears with a probability exponentially small in the number of lattice sites. 

It is interesting to compare the magnitude of the free energy functional for a given 
profile p{x) with that for an equilibrium system that has the same stationary profile. 
This can be realised by coupling a one-dimensional chain of sites along its length to a 
series of particle reservoirs that have their chemical potentials tuned in such a way that 
the desired equilibrium density p{x) in a small interval x G [a, b] is attained. Since the 
equilibrium free energy density is local, the free energy functional for the whole system 
is obtained by summing (6.15) over all these small intervals. This gives 



Now (6.11) reduces to the equilibrium expression (6.21) when a{x) — p{x), the 
equilibrium profile (which is linear). However, the nonequilibrium free energy is obtained 
by maximising the expression (6.11) with respect to a and generally the maximum is 
not realised by ct{x) = p{x). Therefore J-[p] > -^'^^[p] which means that the probability 
of witnessing a particular fluctuation away from the optimal profile for a system driven 
out of equilibrium is suppressed compared to that for an equilibrium system with the 
same optimal profile. This is not always the case however: a similar analysis for the 
asymmetric simple exclusion process (see below) shows that fiuctuations can also be 
enhanced relative to the equilibrium state. 

A final further application of the free energy functional is to explore the optimal 
profiles in the nonequilibrium system after imposing additional global constraints. 
For example, one can ask for the most likely profile given that the overall density 
p — Jq dxp{x) is fixed. It turns out [129] that the the most likely profile has an 



exponential form and, unless the overall imposed density happens to be equal to 
the equilibrium density (in which limit the linear profile is recovered), the density is 
discontinuous at the boundaries x = 0,1. 




(6.20) 




(6.21) 
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6.4- Free energy functional for the partially asymmetric simple exclusion process 
(PASEP) 

Free energy functionals have also been derived for the partially asymmetric exclusion 
process (PASEP) in the case where entry and exit of particles occurs at the left and right 
boundaries respectively, and the bulk bias is to the right (i.e., the forward-bias regime, 
g < 1) [131, 132]. In principle, one could follow through the procedure outlined above in 
Subsection 6.1, where one starts with the full stationary distribution of microstates and 
takes the continuum limit. It turns out to be more straightforward instead to use the 
microscopic distribution to extend the additivity principle discussed in Subsection 6.2. 
Then, by the effective local equilibrium property, one can construct the full free energy 
functional as before. 

When the bulk dynamics are asymmetric, it turns out that two versions of the 
additivity principle come into play. One applies when the left boundary reservoir density 
(t(0) is less than that at the right; the other when the reverse is true. In both cases, 
the additivity relation involves the modified free energy defined by Equation (6.16), but 
where now the current is that found (for example) by an application of the extremal 
current principle discussed in Section 2. That is, 

{min p(l — p) a < a' 
''^["'"'1 , , , . (6.22) 

max p{l- p) a>a' ^ ' 

When (t(0) > o"(l), the additivity formula is the same as for the symmetric case. 
Equation (6.18) (see [132] for details of the calculation). Using the local equihbrium 
property, one then finds the free energy functional to be 

•^m([p]; tx(a), a{h)) = -{b - a) In J(a(0), a{l)) + 

max / dx[p{x)ln{p{x)[l - a{x)]) + {1 - p{x))ln[{l - p{x))a{x)]] . (6.23) 

Again, the companion function a{x) that gives the effective reservoir densities in the 
bulk must match the actual reservoir densities at the boundaries, and must also be a 
nonincreasing function. 

When (t(0) < cr(l), the additivity formula takes on the different form 

^MllMi^^W^f^W) = mill {n[a,y]{[p]](T{a),a{y)) + n[y^b^{[p];a{y),a{b))} ,(6.24) 

Pye{(T{a),a{b)} 

where again we refer the reader to [132] for the details of the calculation. The fact 
that the intermediate reservoir inserted at the point y always takes on a density that is 
equal to that of one of the boundary reservoirs implies that effective reservoir density is 
constant on one side of the point y. Further subdivision on that side then has no effect. 
On the other side the situation is the same as for the previous subdivision: on one side of 
the division, the density will be constant and equal to that of the appropriate boundary 
reservoir, whereas on the other further subdivisions will be required. The upshot of 
this is that the reservoir density function a{x) will be a step function, taking the value 
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(t{x) — a{a) for a < X < y and a{x) — a{h) ior y < x < b for some value y. This value 
will be such that the overall free energy is minimised, that is 

(^(b)) = -(b - a) In J(a(a), a(b)) + 

min < / dx [p{x)ln.{p{x)[l — a{a)]) + {1 — p{x))ln.[(l — p{x))a{a)]] + 

a<y<b { 

rb 

dx [p{x) \n{p{x)[l - a{b)]) + (1 - p{x)) ln[(l - p{x))a{ 

'y 

Note that despite considerable cosmetic differences between this formula and (6.23), 
the two expressions are in fact very similar: in both cases, one needs to find the set of 
reservoir densities (t{x) that leads to an extremum of the same joint functional of p{x) 
and a{x). The fact that in one case, the additivity principle involves a maximum and 
in another a minimum has its origins in the nature of the saddle-point which provides 
the relevant asymptotics [127]. 

We also note that the bias parameter q does not enter explicitly into these equations, 
only implicitly through the relationship between the boundary reservoir densities a{a) 
and cr{b), and the microscopic transition rates a, (3,^,5 and q. Specifically 

a(a) = and a(b) = 1 - , (6.26) 

1 — q 1 — ? 

and so q affects phase boundaries in the a-f3 plane, as was seen in Section 5. The q- 

indcpcndcncc of the free energy functionals docs not contradict the rich structure seen 

in the density profiles in Section 5, since the deviations from the bulk density at the 

boundaries vanish under rcscaling in the N oo limit. This lack of g-dependence also 

implies that (6.11) is not obtained as g — > 1: in other words, the limits N ^ oo and 

g — > 1 do not commute (as we have already seen). The free energy functional in the 

weakly asymmetric limit 1 — q — X/N has been computed explicitly using the matrix 

product approach and on this scale a q (or A) dependence is apparent [133]. 

Closer scrutiny of these free energy functionals for the PASEP [132] shows 
similar properties to those seen for the SSEP above in Subsection 6.3. For example, 
one finds that the stationary profile (here, a constant profile after rescaling in the 
thcrmodynamical limit) is also the most likely profile, except along the boundary 
between the high- and low-density phases along which — as has been discussed before — 
there is a superposition of shocks with the shock location distributed uniformly across 
the system. Then, any one of these shock profiles is found to minimise the free energy 
functional. Again, one can compare the relative size of a fiuctuation away from the 
most likely profile in the nonequilibrium system, and an equilibrium system coupled to 
a reservoir with a spatially- varying chemical potential. In the case where a (a) > (j{b), 
it is found that (as for the SSEP), such fluctuations are suppressed; however when 
(T(a) < a{b) the opposition between the boundary densities and the bulk bias results in 
these fluctuations being enhanced. 

A flnal interesting feature of the PASEP is that density fluctuations in the maximal 
current phase are non-Gaussian. Evidence for this is provided by the existence of a 
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discontinuity at a density of p = | in the probability of witnessing a density p in 
a box located somewhere in the bulk of the system. That the distribution is indeed 
non- Gaussian is confirmed by explicit calculations for the totally asymmetric case 
{q = 0) [132,134] that exploit the relationship to equilibrium surface models outlined in 
Subsection 3.6. 

6.5. Remarks 

In this section we have seen that once one has found a function that is additive when 
two subsystems are connected together via a reservoir with an appropriate density, 
the distribution of density profile macrostates can be constructed given the existence 
of a local equilibrium property. This is an appealing approach, and one hopes that 
it can be applied to a much wider range of systems than exclusion processes. The 
difficulty is, however, that one does not know a priori what form the additivity 
principle shoTild take: here, we have relied on the complete knowledge of the underlying 
distribution of microstates to construct it. Nevertheless, it is worth remarking that 
the free energy functional for the symmetric exclusion process (6.11) has been obtained 
in a purely macroscopic formalism [135]. In this approach it is assumed that, in the 
combined thermodynamic and continuum limit, the macroscopic density profile evolves 
deterministically as a consequence of the law of large numbers. Then, it can be shown 
that the probability of witnessing a deviation p{x) from the most likely profile p{x) in 
the steady state is given by a functional of the most likely trajectory through phase 
space that begins at time t — — oo at stationarity p{x) and is constrained to reach p{x) 
at t = [3]. 

7. Two-species Models with quadratic algebra 

So far we have seen exact matrix product solutions for the steady state of two classes of 
model systems with nonequilibrium dynamics: the open boundary ASEP and PASEP — 
various aspects of which were discussed in detail in Sections 2, 3, 5 and 6 — which had a 
single species of particles hopping on a one-dimensional lattice with open boundaries; a 
two-species model with periodic boundary conditions, which we examined in Section 4. 
In this section we are going to search for 2-species exclusion models that can be solved 
using the matrix product approach. For clarity, we reiterate that in this work the 
number of species n relates to the number of particle species, excluding vacancies. 

It is not known how to perform an exhaustive search of all possible matrix product 
steady states, but it is possible to search a restricted subset where the matrix products 
involved can be systematically reduced using expressions fike (2.39)-(2.41) and (4.4)- 
(4.6). To this end, let us recall the proof of the reduction relations that appfied for 
exclusion models on the ring geometry given in Subsection 4.2. This proof concerned 
stationary weights that were given by a trace of matrices 

/(ti,T2, . . . ,Tiv) = tr(Xr,X,, ■ --X,^) (7.1) 



79 



where the variable Tj = 0, 1, 2, . . . , n indicates the species of particle occupying site i 
{ti = denotes a vacancy) and X^-. is the corresponding matrix whose forms is to be 
determined. We showed that when W{TiTi+i Ti+iTi) is the rate at which particles on 
sites i and i + 1 exchange places, if one can find auxiliary matrices Xi such that 

W{Ti+iTi riTi+i)Xi+iXi — W{Tiri+i ri+iTi)XiXi+i — XiXi+i — XiXi+i (7.2) 

is satisfied, the weights given by (7.1) are stationary. 

If we are to obtain matrix reduction relations, wc require that the auxiliaries X^ 
are scalars (rather than matrices). We must also insist that these reduction relations 
describe an associative algebra, i.e., that no matter what order the reduction relations 
are applied, one always ends up with the same sum of irreducible strings. The various 
ways in which this can be be achieved were formalised, classified and catalogued by 
Isaev, Pyatov and Rittenberg [136]. In this section, we briefly outline the classiflcation 
scheme for the the case of two particle species, n — 2, and show what physical dynamics 
the various possibilities correspond to. We then move on to discuss geometries other 
than the ring. 

7.1. Classification of Isaev, Pyatov and Rittenberg for two- species models 

In this section we will mostly use the notation established above, where vacancies 
are denoted by and particles by 1 and 2. Occasionally, it will be helpful to state 
model dynamics using a natural "charge representation", where particles of species 1 
and vacancies are relabelled as positive and negative charges and particles of species 2 
as vacancies: 

1^+ 2^0 0^-. (7.3) 

Note that our notation differs from that used by [136] . We will also use a shorthand for 
the transition rates 

Wrr' = W{tt' t't) (7.4) 

and, since we are assuming the auxiliary matrices are scalars, we will write them as 

We first restrict ourselves to the case Wi2 > 0, W2q > and Wi2 > 0, i.e. there are 
exchanges in at least one direction between particles and vacancies and between the two 
species of particles. Then, the relations (7.2) become 

WiqXiXq - WqiXqXi = XqXi - XiXq 
W2QX2XQ - WQ2X0X2 = XqX2 - X2XQ 

W12X1X2 - ^21-^2-^1 = 2:2^1 - (7.5) 

It remains to check whether the relations (7.5) arc consistent. The approach of 
Isaev, Pyatov and Rittenberg [136] (see also that of Karimipour [137]) is to generalise 
the reduction process leading to (2.45). That is, one seeks to use (7.5) to reduce any 
product of matrices to a sum of irreducible strings. First we fix the order of irreducible 
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strings as XqX^X". Then we require that reducing an arbitrary product U of the 
matrices Xi to a sum of irreducible strings 

U^Yl ai,m,nX!,X^X^ (7.6) 

l,m,n 

is independent of the order in which the reduction is carried out using the elementary 
rules (7.5). Reduction rules which satisfy this requirement are referred to by Isaev, 
Pyatov and Rittenberg as PBW-type algebras [138]. 

For the two-species case the requirement is that the two possible reductions 
of the product X1X2XQ, illustrated schematically below, give equivalent results. 

X1X0X2 XqXiX2 . 



XIX2XQ ^ X0X2XI 

X2XIX0 — ^ X2X0XI 
Calculating explicitly using (7.5) yields [136] 

Xi Wo2{wi2 - W21 - Wio + Woi)^0^2 + X2 Woi{wi2 - W21 + W20 - ^02)^0^1 

+ xo W2i{w2o - W02 - wio + i(;oi)^2-'^i + a;ia;2(tyi2 - W21 + W20 - wio) Xo 
+ a;iXo(u'2i - W02) X2 + X2Xo(wio - wu - W20 + W02) Xi ^ , (7.7) 

which results in the following six conditions on the hopping rates Wij from the 
requirement that each of the six terms in the above equation vanish: 

XiWo2{Wi2-W2l-Wio + Woi) =0 , (7.8) 

X2 W0l(?«12 - W21 + W20 - Wq2) = , (7.9) 

XoW2l{w20-Wo2-Wio+Woi)=0, (7.10) 

XiX2{wi2 - W21 + W20 - Wio) = , (7-11) 

a:;ia:;o(w2i - W02) = , (7-12) 

2;22;o(«'lO - Wi2 - W20 + W02) = . (7.13) 

The solutions of these equations were classified in [136]. Here we summarise the various 
nontrivial solutions which are physically relevant and the corresponding models which 
have been previously studied in the literature. The classification hinges on how many 
of the Xi are zero. 

7.1.1. Class A (includes Karimipour's model of overtaking) The first class of solutions 
has none of the Xi = 0. There are then two possible solutions to (7.8)-(7.13). 

• Solution AI Wij = w Physically, this corresponds to symmetric exclusion 

with two species which, although labelled distinctly, have identical dynamics. 
Matrices are easily found by setting X2 = {x2/xi)Xi which reduces (7.5) to the 
single species condition w{XiXq — XqXi) = XqXi — XiXq. On the ring this system 
has a simple steady state where all allowed configurations are equally likely. 
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• Solution All Wij = Vi — Vj > for i < j and Wij — ior i > j. This corresponds 
to a model where the hop rates are totally asymmetric and of the form 

10^01 20^02 12"^-^"^ 21. (7.14) 

That is, since Vi > V2, a. particle of species 1 is faster than a particle of species 2 
and overtakes with rate Vi — V2- The corresponding matrix algebra is 

ViXiXq — XqXi — XiXq 

f 2X2X0 = 0:0X2 - X2X0 

{vi - V2)X,X2 = X2X1 - X1X2 . (7.15) 

This solution of (7.8)-(7.13) was first noted and studied by Karimipour [137] and 
followed up in [139]. As it happens, this model can be generalised to more than 
two species, and as such will be discussed in its full generality in Section 8.2. 
Representations are given in Appendix D Equation (D.19). 

7.1.2. Class B (includes asymmetric second-class particle dynamics) This second class 
of solutions has one of = 0. The first two famihes of solutions within this class are 
obtained by taking X2 — 0, and then choosing the hop rates to satisfy (7.8), (7.10) and 
(7.12) in two different ways. 

• Solution BI W12 — W20 — P2, 1^21 — W02 — q2 , Wio — pi and Wqi — qi with 
P2 — q2 — Pi — qi- This corresponds to a model with the hop rates 

pi P2 P2 , ^ 

10^01 20^02 12^21 (7.16) 

91 92 ?2 

or, in the charge representation, 

+ -^-+ 0-^-0 +0^0+ . (7.17) 

qi <?2 92 

This model is thus an asymmetric generalisation of the second-class particle of 
Section 4 and was first used by [78] for the special case q2 = qi, P2 = Pi- The 
corresponding matrix algebra is 

PiXiXo - giXoXi = xqXi - xiXq 
P2X2X0 — ^2X0X2 = X0X2 

P2X1X2 - 92X2X1 = - X1X2 (7.18) 
A representation of this algebra is given in Appendix D. 

• Solution BII W02 = W21 = 0, W12 — ck, Wiq — p, Wqi — q and W12 — /3. The 
corresponding model has hop rates 

10;^ 01 20A02 12^21 (7.19) 

g 

or, in the charge representation, 

+ -^-+ 0-A-O +o4o+ . (7.20) 
9 
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This model is another generahsation of the second-class problem first studied with 
p = 1, g = for a single species two particle [90] and later generalised in [107]. The 
matrix algebra is 

pXiXo - qXoXi = xqXi - xiXq 

aX2Xo — X0X2 

pXiX2 = - X1X2 . (7.21) 

This algebra can be mapped onto that used to solve the PASEP with open 
boundaries, (5.1)-(5.3), if one takes 

X2^\V){W\ (7.22) 

and Xq — 1, Xi — —1. Thus one can make use of results from Section 5, along with 
techniques for models on the ring described in Section 4 to analyse various cases of 
this model [107-112]. 

Choosing xi or xs to be zero results in the same solution as BIT under relabelling 
particles, and one additional solution when xi — 0. This is 

• Solution Bill woi — W21 — 0, Wio — a, W21 — /3, ^20 = P, W02 — q and 

p — q — a — p. The corresponding model has hop rates 

lOAOl 20;^02 12^21 (7.23) 

and the matrix algebra is 

aXiXo = xqXi 

PX2X0 - qXoX2 = X0X2 - X2X0 

pX,X2 = X2X, . (7.24) 

However on the ring this results in a steady state in which all allowed configurations 
arc equally likely: as can be seen from (7.24) in any periodic string of matrices, all 
X2 and Xq can be eliminated through the first and third relations. 

7.1.3. Class C (includes non- overtaking two-species dynamics) The third class of 
solutions is obtained by taking two of the scalar quantities Xi equal to zero. Equations 
(7.11)-(7.13) and two of (7.8)-(7.10) are then automatically satisfied. The equation that 
remains has a structure that is independent of which of the Xi arc taken to be nonzero, 
so we take Xi — X2 — Q which leaves us to satisfy (7.10). This can be done in two ways. 

• Solution CI wiQ = pi, wqi = gi, W20 = P2, W02 = ?2 and pi - qi ^ P2 — q2- This 
corresponds to a model in which all six exchanges may take place 

Pl P2 W12 

10^01 20^02 12^21 (7.25) 

91 92 W21 

and that has the matrix algebra 

PiXiXo - qiXoXi = XqXi 
P2X2X0 — g2^o^2 = X0X2 

W12X1X2 - W21X2X1 = . (7.26) 
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Unfortunately this algel)ra is not useful in describing a physical system with periodic 
boundary conditions. To see this one can take for example a representation where 
Xo = 1 and xq = Pi — Qi = P2 — (I2 which satisfies the first two relations of 
(7.26). However the third relation of (7.26) has the form of a vanishing deformed 
commutator 

X1X2 - rX2Xi = . (7.27) 

Thus if we commute say a matrix Xi all the way around the ring we will end up 
with the same matrix product multiplied by a factor r^^ where N2 is the number 
of species 2 particles on the ring and r = W2i/wi2- Thus, only for W21 = W12 
(r = 1) can we use this algebra on a periodic system; alternatively, it can be used 
for general r on a closed segment as we discuss in Section 7.2.2. 

• Solution CII Equation (7.10) can also be satisfied by taking wu = W21 = and 
unrestricted choices for the remaining rates Wio = Pi, Wqi = qi, W20 = P2, "^02 = 52- 
This has similar dynamics to the previous model, but with no exchange of species 
1 and 2 

Pl P2 

10^01 20^02 (7.28) 

The matrix algebra is similar to (7.26) but with the third relation absent. In 
contrast to the previous case, this algebra can be used to describe the above 
dynamics on a ring [140]. This model and its generalisation to multiple species 
will be discussed further in Section 8. 



7.1.4. Class D (includes the cyclically symmetric ABC model) In this case, we take all 
Xi — Q which leaves all six hop rates free. 

WIO "'20 W12 ,_ 

IOf^OI 20^02 12^21 (7.29) 

WOl "'02 W21 

The algebra is a set of deformed commutators 

w^oXiX^ - WoiXoXi = (7.30) 

t(;20-'^2-'^0 - W02-'^0-'^2 = (7-31) 

^/;i2XiX2 - W21X2X1 = . (7.32) 

Explicit forms for Xq, Xi and X2 are discussed in Section 7.2.2. On a periodic system, 
as with case CI discussed above, we require that a matrix product is left invariant after 
commuting one of the matrices once around the ring. The condition for this is that 



'woi 


no 


W2I 


n2 


W02 


no 


W12 


ni 




ni 


W20 






_Wi2. 




_W20_ 




_W21_ 




Wol. 




_Wo2_ 



n2 



(7.33) 



where rii is the number of particles of species i. These conditions are satisfied for 
example by the ABC model [141,142] in the special case where all particle numbers are 
equal Uq = rii = n2. The fact that the rhs of (7.32) are all zero implies that to use 
these algebraic relations, detailed balance must hold in the steady state. Therefore the 
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condition (7.33) is actually the condition for detailed balance to hold. The corresponding 
energy turns out to be an interesting long-range function [142] . 

In the ABC model, the particles and vacancies are relabelled as 

l^A 2^B O^C (7.34) 

and have dynamics 

AB^BA BC^CB C A^ AC (7.35) 

111 ^ ' 

where we take q <1. This corresponds to the choice of rates Wi2 = "^^20 = "^oi = Q and 
W21 = = WiQ = 1. Note that the dynamics are invariant under cyclic permutation 
of the three particle types, and that A particles tend to move to the left of B particles, 
B particles tend to move to the left of C particles and C particles tend to move to 
the left of A particles. In the steady state of this model this results in a strong phase 
separation into three domains of A, B and C in the order ABC. In the limit N ^ 00 
with g < 1 fixed, the domains are pure in that far away from the boundary of the 
domain the probability of finding a particle of a different species (to that of the domain) 
tends to zero. That is, particles from a neighbouring domain may penetrate only a finite 
distance. In the case where we have equal numbers of each species condition (7.33) is 
satisfied and one can use the matrix product to calculate the steady state exactly [142]. 
This model also exhibits an interesting phase transition in the weakly asymmetric limit 
where q varies with system size N as q — exp{—l3/N), thus approaching unity in the 
thermodynamic limit. Then, according to the value of /3, the steady state can either 
order into three domains which are rich in A B and C (but are not pure domains) or a 
disordered phase where the particles are typically in a disordered configuration [143] . 

7.2. Geometries other than the ring 

In the previous subsection, we summarised a classification of all the possible sets of 
dynamics for which a matrix product state with scalar auxiliaries (i.e., one that has a 
set of reduction relations) exists with a unique decomposition into irreducible strings 
of matrices. We also assessed whether these matrix product states could be used for 
models on the ring. We now extend this enquiry to other one-dimensional geometries 
that can be constructed. 

7.2.1. Open boundary conditions We first consider open boundary conditions, i.e., 
those where particles can enter and leave at the boundaries. The most general dynamics 
is to have rates aij at which a particle of species i (or a vacancy if i = 0) is transformed 
into a particle of species j at site 1 and rates at which a particle of species i is 
transformed into a particle of species j (or a vacancy of j — 0) at site N. If we are to 
have statistical weights of the form (2.35), i.e., 

fin, . . . ,T;v) = {W\X,,Xr, ■ ■■Xr.lV) (7.36) 
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we then obtain additional conditions involving the matrices and vectors {W\ and \V). 
These may be written as 



where an = - aij and A, = - (3ij. Since Y,j o^ij = Aj = we require 



In the open boundary case, typically, in addition to the conditions for the solution 
classes of Section 7.1, one has further constraints on the boundary rates to allow 
conditions (7.37,7.38) to be satisfied. We do not attempt to catalogue all solutions 
here, rather we point out which of the solution classes of Section 7.1 may, in principle, 
be used in an open system and reference some examples. Some general considerations 
of the open boundary conditions for which one has a matrix product state using the 
quadratic algebra (7.5) are discussed in the two species cases in [144]. 

• Class A Solution AI, which corresponds to symmetric exclusion of two differently 
labelled species with identical bulk dynamics, can be used in the open boundary 
case where the two species are injected and extracted with different rates, under 
certain conditions. Solution All can be combined with certain open boundary 
conditions that were found by Karimipour [139]. We shall specify these boundary 
interactions and discuss this model more thoroughly in Section 8. 

• Class B Both solutions Bl and Bll can be used in models with open boundary 
conditions . One case that has been studied is the dynamics implied by solution Bll 
with q = [18, 145-147]. This is sometimes referred to as the bridge model because 
one has two particle species with opposite velocities that slow down to exchange 
places when they meet, rather like cars on a narrow bridge. Unfortunately, the 
full phase diagram of the bridge model — and in particular an interesting broken 
symmetry region — cannot be described by a quadratic algebra (i.e., one where 
the auxiliaries are scalars and reduction relations are implied on the steady-state 
weights). The BII case has also been studied [148, 149]. 

• Class C Recall that Class C comprises solutions where only one of the Xj are 
nonzero. Hence, it is not possible to satisfy that sum rule (7.39), and hence the 
dynamics implied by these solutions cannot be solved in an open system using a 
simple matrix product algebra. 

• Class D Algebras within class D can only be used in models with open boundaries 
under the special conditions of detailed balance, which generally does not hold. On 
the other hand, as we discuss below, detailed balance does apply in conservative 
models with closed boundaries and thus such algebras are relevant there — see below. 




(7.37) 




(7.38) 




(7.39) 
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Figure 29. Typical configuration in the steady state of the partially asymmetric 
exclusion process on a closed segment of N sites with q < 1. 



7.2.2. Closed segment A closed segment is a one-dimensional lattice of size N with 
dosed boundary conditions. That is, particles cannot enter or leave at the left boundary 
or right boundaries (sites 1 and A^). For the moment, let us consider conserving models 
where particles are neither created nor destroyed in the bulk: the only moves that are 
allowed are those where particles or vacancies on neighbouring sites exchange places. It 
is straightforward to show using a Kolmogorov criterion (6.2) that detailed balance is 
satisfied in the steady state of any such model. The reason for this is that in order to 
create a loop in configuration space, every exchange of a pair must be accompanied by an 
exchange in the opposite direction (as long as the dynamics is not totally asymmetric). 
Hence the forward loop contains the same set of exchanges as the reverse, and (6.2) is 
satisfied. Some quadratic algebras for closed segments and segments closed at one end 
and open at the other were discussed in [144]. 

The simplest example of this type is the partially asymmetric exclusion process 
(PASEP) on the closed segment i.e. just one species of particle [150]. The appropriate 
matrix algebra in this case is 



Taking g < 1 the configuration with the highest weight will be the configuration with all 
particles stacked up to the right end of the lattice as shown in Figure 29. The weight of 
any configuration can then be obtained by moving particles to the left from the maximal 
weight configuration and multiplying the weight by a factor q at each move as implied 
by (7.40). Then the statistical weight of any configuration ti,T2, . . . ,tn will be given by 



The steady state weight (7.41) is in fact a Boltzmann weight with an energy function 
H oc I lng| X]fc ^''"fc the dynamics respects detailed balance with respect to this 
weight. 

Let us now turn to Class D models, where the matrix algebra takes the form 
of three deformed commutators (7.32). An explicit form for the operators Xi,X2,X3 
involves tensor products of matrices which obey deformed commutation relations 



DE - qED = . 



(7.40) 




(7.41) 



and the maximal weight is 




Xi = E{l02i/uJi2) ® 1 

X2 = \®D(S> E{UJ02/UJ20) 

Xo = ^(cjio/c^oi) ® 1 ® 



(7.42) 
(7.43) 
(7.44) 
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where 

DE(r) - rE{r)D = . 
So, for example, 

[wioDE{uJio/uJoi) - woiE{uJio/uJoi)D] E{uj2i/uji2) 



(7.45) 



(7.46) 



To obtain statistical weights here, we need to prescribe a suitable contraction 
operation i.e. appropriate boundary vectors. Using the representation of the deformed 
commutator algebra in Appendix D (D.29) we see that if there are ni species 1 particles, 
n2 species 2 particles and Uq vacancies then 



{w\ = (o| (g) (o| (o| \v) = Im) 



I no) 



(7.47) 



will give a nonzero contraction for those configurations with the correct number of each 
species. In particular for the reference configuration the weight is 

— ] . (7.48) 

7.2.3. A single defect on an infinite system In order to investigate the structure of 
shocks, Derrida, Lebowitz and Speer studied a single second-class particle on an infinite 
system [151]. The idea was to describe the stationary measure as seen from the second- 
class particle by a matrix product. The model considered was the partially asymmetric 
generalisation with dynamics 



10^01 

q 



20^ 02 

q 



12^ 21 

q 



(7.49) 



Taking the second-class particle to be at the origin and considering a window of m 
sites to the left and n sites to the right of the second-class the stationary probabilities 
(note, not weights) are written as 

-1 



P{r-m, Ti, . . . , Tn) = {W\ 



A 



1^) 



(7.50) 



where A is once more the matrix corresponding to a second-class particle. In order for 
this form to hold for a window of arbitrary size i.e. all values of n, m = . . . oo one 
requires that 

{W\C^{W\ 

C\V) = \V) (7.51) 

where C — D + E. To ensure that the probabilities are correctly normalised we further 
require 

{W\A\V)=1. (7.52) 

The bulk dynamics falls within the solution class BI, and so the matrix algebra 
(7.18) should be used with Pi — P2 — P and Qi — q2 = P- This particular application has 



88 



the unusual property that the choice of Xq and Xi has physical consequences. This is 
because Xq and Xi cannot be scaled out of equations (7.51). The correct choice is made 
by insisting that the desired asymptotic densities p+ far to the right of the second-class 
particle (n — > oo) and p_ far to the left (m — > oo) are obtained. It turns out that one 
should take 



Choosing p+ > p_ corresponds to a shock profile as soon from the second-class particle 
with the second-class particle tracking the position of the shock. The structure of the 
shock was analysed, in particular the decay to the asymptotic values is exponential. An 
interesting result is that the characteristic length of the decay becomes independent of 
the asymmetry when 



To obtain a representation of the required matrices and vectors, one can use (D.21)- 
(D.23) given in Appendix D for solution class BI. Then, using (7.51), one can construct 
(1^1, \V) to be eigenvectors D + E with eigenvalue 1. The expressions for the boundary 
vectors are, however, quite complicated in this representation. An alternative approach, 
used in [151], exploits an infinite-dimensional (rather than semi- infinite) representation. 
Finite dimensional representations along special curves were studied in [152]. 

7.3. N on- conserving two-species models with quadratic algebra 

In the foregoing, we have assumed that the bulk dynamics conserve particle numbers. In 
this section we show how the matrix algebras which were encountered can be adapted to 
models that have non-conservative dynamics. The basic idea is to augment the dynamics 
with additional processes that move between different sectors (i.e., configurations with 
a particular number of particles) . These moves will be constructed in such a way that 
detailed balance between sectors is satisfied, even though detailed balance does not hold 
within a sector. As such, one can realise any prescribed distribution of sectors: we 
choose weights that correspond to the grand canonical ensemble that was introduced in 
Section 4 as a means to study systems with fixed particle numbers. We first discuss this 
method in detail, and then give two concrete examples. 

7.3.1. N on- conserving dynamics which generate a grand canonical ensemble As we 
have just described, the aim is to construct some dynamics such that the statistical 
weight of a configuration with M particles (one a designated species) on an A^-site ring 
is 



where are the matrices that give the steady state for some conservative process, and 
M is a fugacity that controls the mean particle number M. Let us denote transition 



Xo = {p-q){l- p-){l- p+) 
xi = -{p-q)p-p+ 



(7.53) 
(7.54) 




(7.55) 




(7.56) 
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rates for this conservative process as Wo{C C), and the additional rates that serve 
to increase or decrease M by one as W±{C C). The statistical weights must then 
satisfy the master equation 

[/(C; M)W,{C' ^ C) - /(C; M)W^{C ^ C')\ + 



[/(C; M - l)W+{C' ^ C) - /(C; M)W_{C ^ C')] + 

c 

J2 [fiC; M + l)W.{C' ^ C) - /(C; M)W^+(C ^ C')] = . (7.57) 
c 

Now, the first summation is what appears in the master equation for the conservative 
process, and since all the weights for fixed M are proportional to , this first 
summation vanishes. The second and third summations can be made to vanish if we 
choose 

W4C ^ C) _ f{C'- M - 1) _ 1 tr(X,.X,. ■■■^r'^. 



W^{C'^C) f{C;M) utT{Xr,Xr,---Xr^) ^^-^^^ 

where Tj and r/ are the site occupation variables specifying to configurations C and C 
that have M and M — 1 particles respectively. As we can see, this takes the form of a 
detailed balance relation (6.3). The two examples that follow serve to demonstrate. 

7.3.2. Grand- canonical dynamics for the ASEP with a second-class particle Recall that 
in Section 4 we considered a version of the ASEP with a second-class particle, i.e., a 
model with dynamics 

10-i>01 20^02 12^21 . (7.59) 

This falls within the family of solutions denoted BII above. In the analysis of this model 
we introduced a fugacity m as a trick to expediate calculation of statistical weights, 
tuning its value to yield the desired mean particle density on the ring in the canonical 
ensemble. 

This ensemble can also be generated using physical dynamics as we have just 
described. This is achieved by noting that the algebra Let M count the number of 
positive charges (species 1 particles). We note that the algebra (7.21) imphes 

iTj-.-X^X,...) _ (3x0 . . 

tr(---A:iA:2---) axi' ^ ' ' 

Hence from (7.58), we see that the grand canonical ensemble is realised physically if we 
introduce two non-conserving processes 

+ oio - (7.61) 

7 

where the rates satisfy 

^ = -^^-^ . (7.62) 
7 au xi 
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We remark that S and 7 as free parameters that enforce a particular fugacity u in 
the grand canonical ensemble generated by these non-conservative dynamics. This 
particular implementation of a model with non-conserving dynamics was used in [153] 
to generate a grand canonical ensemble dynamically. 

7.3.3. A model in which no particle numbers are conserved A second non-conserving 
generahsation of the ASEP on a ring with a second-class particle is best described using 
the charge representation. For simplicity, we take the dynamics of BII with p — 1 and 
q — and take a — P — 1. The case q was studied in [154]. Thus, the conserving 
dynamics are 

+ 0^0+ 0-^-0 +-^-+ . (7.63) 

To these dynamics wc seek to add the following non- conserving moves 

+ 0;^00 0-;^00. (7.64) 
11 

in such a way that the BII matrix algebra (7.21) can be exploited. To achieve this, we 
require matrices such that (7.58) is satisfied for the non-conserving moves. That is, 

l tr(---XoXo---) _ l tr(---XoXo---) 

utT{---X+Xo---) Mtr(---XoX_---) ^ ■ ' 

where it is a fugacity counting the total number of positive and negative charges. To 
make contact with more familiar problems, let us put — D, X_ — E and Xq = A. 
We see that the previous equation is satisfied if 

u = - and A^^A, (7.66) 

w 

i.e., we require A to be a projector. This is achieved using the matrices for the second- 
class particle problem of Section 4, (4.4-4.7) i.e. A = where D\V) = \V), 
{W\E — {W\ and = 1. The statistical weight of a configuration is then given by 

fin, r2, . . . , tn) = tr{Xr,X^, ■ ■ ■ r^v) (7.67) 

where M is the number of vacancies in the configuration. 

In [154] it was shown that an interesting phase transition arises as w is varied. For 
w < 1 there is a tendency to create positive and negative particles and this ultimately 
leads to a vanishing steady-state density of vacancies. On the other hand for w > 1 
there is a tendency to eliminate positive and negative charges at the left and right 
boundaries respectively of domains of vacancies. Thus domains of vacancies tend to 
grow and in the steady state this results in a finite density of vacancies. The transition 
is remarkable in that it is a phase transition in a periodic one dimensional system with 
local dynamics without an absorbing state. Here, we shall use the generating function 
technique of Section 3.2.3 to quickly obtain the asymptotics of the partition function 
and hence demonstrate the phase transition. 
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First, notice from (7.64) that configurations with no vacancies are dynamically 
inaccessible. Therefore all configurations have at least one vacancy and we write the 
partition function as 

Zn = tr [wA(D + E + wA)^'^] = w{W\ (D + E + wA)^-^\V) (7.68) 

(An unimportant subtlety is that we have ignored the degeneracy factor in placing a 
given configuration of particles on the periodic lattice; this factor is bounded from above 
by N and hence does not contribute to the exponential part of Zjsf which controls the 
macroscopic physics, as will be seen below.) 
We introduce the generating function 

oo ^ 

F{z) = J2 ^""Z^ = MW\^^-^\V) (7.69) 

N=l 

where we define G — D + E + wA. Now, we may write 

1 - zG ^ (1 - rjD - uA){l - rjE - uA) (7.70) 

where 

z — 77(1 — 7]) zw = 2v{l — Tj) — . (7-71) 
Then using the generating function technique and (7.71) we obtain 

F(z)^wz{W\- ^ ^ -\V) (7.72) 

- , = 1 ^1 • (7-73) 

(1 — ?7 — vy (1 — -^Y — wz 

The singularities of F{z) are a square root singularity at z — 1/4 (coming from 
ri{z) = (1 — Vl — 42;)/2) and a pole at z — w/ {l + wY which only exists if w > 1. Using 
the results of Appendix B one quickly obtains the asymptotic behaviour of 

1 w 

a w<l ^iv~^^3^^^ (7.74) 



w-l{w + If" 
w + 1 w 

Due to the form (7.67) of the weights, 9, the density of vacancies, is given by 



it u,>l Ziv ~ . . 7 : ' ■ (7.76) 



^ lim = lim , (7.76) 

where M is the average number of vacancies in the system. Therefore we obtain 

e = for w<l (7.77) 
w — 1 

9 = for w>l (7.78) 

w + 1 ^ ' 

demonstrating the phase transition in the density of vacancies at w — 1. 
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8. Multispecies Models with quadratic algebra 



In the previous section we identified all possible two-species dynamics that have a steady 
state of matrix product form with a quadratic algebra i.e. obeying relations of the type 
given by (7.5). We deferred the discussion of two families of dynamics, as both can 
be generalised in a straightforward way to arbitrarily many species of particles. We 
now examine these cases in more detail. Further discussions of multi-species quadratic 
algebras are given in [136,155-158]. 

8.1. ASEP with disordered hopping rates 

We first consider a variant of the ASEP on the ring in which each particle has a 
different forward and reverse hop rate (and there is no overtaking) [140,159]. This is a 
generalisation of the two species model, case CII of Section 7. There are M particles 
and N — M vacancies on the ring of size N . We give each of the M particles an index 
/X, and introduce for each fi a pair of rates and so that the dynamics of particle 
are 

/xO^O^i. (8.1) 

The corresponding matrix algebra is a generalisation of class CII of Section 7 to M 
particle species. It takes the form 

p^D^E - q^ED^ = (8.2) 

where we have adopted the notation Xq = E and = D^. One possible representation 
of these matrices is given in Appendix D as Equation (D.14). However, it is simplest to 
perform calculations not with a representation, but directly from the algebra (8.2). 

On a periodic system let the steady state weight of a configuration specified by 
{riju} where is the number of empty sites in front of particle be fnini, . . . jTIm)- 
The matrix product form for /jv is 

/Ar(ni, 712, riM) = Tr [D^ E^' E^' . . . Dm E^^] (8.3) 

Using (8.2) ioi n — M m (8.3) gives the relation 

Q.M 1 

/jv(ni, . . . , um) /iv(ni, . . . , riM-i + 1, njw - 1) = — fN-\{ni, nu-i, um - !)■ 

Pm Pm 

The procedure is continued using, in sequence, (8.2) with /i — M — 1, M — 2 . . .1 to 
obtain 



M-l 
^_ J-J 







/jv(ni, . . . ,nM) 



'M-l . M 



fN-i{ni, . . . ,nM - 1) 



The effect of this manipulation has been to commute a hole initially in front of the 
particle M backwards one full turn around the ring. The result is that the weight of a 
configuration of size N is expressed as a multiple of the weight of a configuration of size 
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N — 1 with one hole fewer in front of particle M. Repeating the commutation procedure 
for a hole initially in front of particles M — 1, M — 2, ... 1 implies that the weights are 



M 



/iv({n^}) = JJ^;:Vm(0,...,0) where 
n=i 

and we can take /m(0, . . . , 0) = 1. 



'M~l fJ. 

n i 



j=^i+l-^ 



M 

n- 

k=l 



(8.4) 



8.1.1. Bose-Einstein Condensation An interesting phenomenon that may occur in 
this system is that of condensation. In that case the steady state is dominated by 
configurations where one particle has an extensive value for [140, 160]. This is easiest 
to understand when = i.e. we consider only forward hops. Then (8.4) reduces to 
= 1/Pfj, and we may write 

fNiirif,}) = exp I - ^ n^e^ J where = Inp^ . (8.5) 

The weight (8.5) clearly has the form of the weight of an ideal Bose gas with kT — 1. 
Here the 'bosons' are vacancies and the Bose states correspond to particles with the 
energy of the state determined by the hop rate of the particle. The equivalent of the 
density of states for the Bose system is the distribution of particle hop rates which we 
denote cr(p). If there is a minimum hop rate Pmin this will correspond to the ground 
state of the Bose system. Then, if for low hop rates the distribution of particle hop 
rates follows 

a{p) ~ (p - Pminy with 7 > , (8.6) 

Bose condensation will occur for high enough vacancy density or equivalently low enough 
particle density. When condensation occurs the slowest particle (the one with the 
minimum hop rate Pmin) has an extensive number of vacancies in front of it whereas 
the rest of the particles will have gaps in front of them comprising some finite number 
of vacancies. Thus a 'traffic jam' forms behind the slowest particle. In [12] it is shown 
that the ASEP on a periodic system may be mapped onto a Zero Range Process and 
the condensation transition is fully discussed in the context of the Zero Range Process. 

Recently, further progress has been made in understanding the case where Q/j, ^ 
by using extreme values statistics and renormalisation arguments [161]. 

8.2. Karimipour's model of overtaking dynamics 

This is the generalisation of the two-species case All, which had the two particle species 
hopping with different speeds, and the faster ones overtaking the slower. The many- 
species generalisation allows an arbitrary number of particle species labelled n = 1 . . . P. 
A particle of species /i has a forward hop rate or 'velocity' v^, and adjacent particles 
and 1/ exchange places with rate V/j, — Vi, if V/j, > v^,, i.e., 

/xO^ 0/x 

fiv ^-^^ V H if v^> . (8.7) 
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The model provides an elegant generalisation of the open boundary TASEP and its 
algebra. The generahsation of (7.15) to the multispecies case is 

V/^Xi^Xq = xqX^ - Xf^Xo II > 

- v,,)XnX^ = Xi,X^ - x^j^X^ Vi^> v„ . (8.8) 

Writing Xq — E, X^ = D^/P and making the convenient choice xq — 1, x^ — —v^/P 
for // > yields 

D^E = —D^ + E 

D^D,^ —^D, Vi.>Vi^- (8.9) 

On a periodic lattice the relations (8.9) are the only ones which need be satisfied and 
this can be done by choosing a scalar representation E = 1/e and = where e 
is chosen to be less than the lowest hop rate. Thus on a ring all allowed configurations 
are equally likely. 

More interestingly, the relations (8.2) may be used for an open system with suitable 
boundary conditions. As discussed in Section 7.2.1 we require X];^i=o'^m ~ ^ which with 
the above choice of xo and implies 

At the left hand boundary a particle of species is injected with rate cc^u if the first site 
is empty. The choice — av^/ P reduces (7.37) to 

{W\E^-{W\ (8.11) 
a 

At the left hand boundary a particle of species /i leaves the system with rate (5^. The 
choice = + P — 1 ensures that explicit representations of the matrices can be 
constructed (see Appendix D) and (7.38) becomes 

Relations (8.9), (8.12) and (8.11) provide an elegant generalisation of the usual ASEP 
relations (2.39)-(2.41) which are recovered when we have just one species of particle with 
hop rate vi — 1. In general, an infinite dimensional representation of (8.9) is needed to 
satisfy the boundary conditions (8.12); see, for example, (D.19) in Appendix D. 

The fully disordered system is realised when each particle that enters has a velocity 
drawn from some distribution cr(v), with support [vmiiDOo] so that the lowest velocity 
is t'min, and (8.10) is replaced by 



dva{v)v^l. (8.13) 



Thus, when site one is vacant a particle enters with rate a and its velocity v is assigned 
according to the distribution ij{v)v giving an effective rate a{v) = (T{v)va. Therefore 
we can think of a reservoir of particles with velocities distributed according to a{v) 
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each attempting to enter with rate va. A particle reaching site N exits with rate 
(3{v) = V + P — 1. Then we can replace in the discrete case by X[v) in (8.8). Letting 
Xq = E and X{v) = a{v)D{v), xq = 1 and x{v) = —a{v)v recovers the same algebra 
(8.9), (8.12) and (8.11) with replaced by D{v), e.g., 

D(v)D(v') = -^^D(v') —D(v) for v > v' . (8.14) 

V — V V — V 

Karimipour [139,162] studied, amongst other things, the phase diagram of the open 
boundary model and found the structure depends on the distribution a{v) through the 
parameter 

which characterises the form of a{v) near the lowest hop rate I'min- If ^[o"] < the three 
phases (low-density, high-density and maximal current) of the open boundary ASEP 
remain although the phase boundaries depend on parameters such as l[a] and I'min- 
However if l[a] > the high-density phase is suppressed and only the low-density and 
maximal current phases exist. 



9. More complicated matrix product states 

Recall that in Section 3.1.2 we set out a general cancellation scheme for matrices that 
guarantees a stationary solution of a master equation for a process on the ring. This 
involved a set of ordinary matrices Xt and their auxiliaries Xr- By restricting to the 
case of scalar auxiliaries, we showed, in the case of two particle species in Section 7, 
that it was possible to determine all possible sets of dynamics that give rise to a matrix 
product steady state. In this section, we shall consider the more complicated case where 
the auxiliaries are matrices or more complicated operators. 

Here it is not possible — as far as we are aware — to perform an exhaustive search 
for solutions. However, it has been shown that if one has a one-dimensional lattice with 
n particle species, open boundary conditions and arbitrary nearest-neighbour dynamics 
in the bulk, there does exist a matrix product solution, involving auxiliaries that in 
general will not be scalars, obeying the algebraic relations set out in Section 3.1.2. In 
the next subsection we review this existence proof which is due to Krebs and Sandow [92]. 
Unfortunately, this proof does not lead to any convenient reduction relations, like (2.39)- 
(2.41) for the ASEP, that might be in operation. Also, the proof is not constructive in 
that it requires the steady state to already be known in order to construct explicit 
matrices. Rather, the proof demonstrates that there are no internal inconsistencies in 
the cancellation mechanism of Section 3.1.2. 

In the last two subsections we discuss two models of which we are aware that give 
examples of a matrix product state with matrix auxiliaries. 
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9.1. Existence of a matrix product solution for models with open boundaries 

The existence proof of a matrix product state [92] (see also [163] for the discrete time 
case) applies to models with open boundary conditions and arbitrary nearest-neighbour 
interactions. Let us restate the master equation for a general process in the form 
introduced in Section 3.1.2. It reads 

^/(n, . . . , Tjv) = Hf{Ti, ...,tn)= ( + + hRj fin, ...,tn), (9.1) 

where the operator hi^i+i applied to the function of N state variables Tj generates the 
gain and loss terms arising from interactions between neighbouring pairs of sites in the 
bulk, and fiL and fiR do the same at the boundaries. Here we shall consider models that 
have n distinct particle species in addition to vacancies: i.e., the occupation variables 
take the values Ti — 0,1, ... ,n. The matrix product expressions 

fin, T2, • • • , tjv) = {W\Xr,Xr, ■ ■ ■ X^JV) (9.2) 

provide the steady state solution to (9.1) for some set of matrices Xr and vectors {W\ 
and \V) if one can find a second set of auxiliary matrices X-r such that the relations 

k,+,{W\ ■ ■ -X^^Xr^^, ...\V) = {W\--- [Xr.X^^^, - X^X.J ■■■\V) (9.3) 
hL{W\Xr, ■■■\V) = -{W\Xr, ■■■\V) (9.4) 
hn{W\---X,jV) = {W\---X,jV) (9.5) 

hold. Then, one obtains a zero right-hand side of the master equation (9.1) via a pairwise 
cancellation of terms coming from (9.2)-(9.5). Thus, relations (3.21-3.23) for the ASEP 
generalise to 

hi,i+lXriXri_^_i = Xt--Xt-._^^ - Xr^Xr^_^_J^ (9.6) 
hL{W\Xr, = -{W\Xr, (9.7) 

hnX^^lV) = Xr^\V) . (9.8) 

Relations (9.6-9.7) give the general form for a cancellation mechanism involving matrix 
(or possibly tensor) auxiliaries. Such a cancellation scheme has been proposed in 
various contexts [26, 164, 165] including a generalisation to longer (but finite) range 
interactions [166]. Some consequences of this scheme have been explored in [167]. 

The existence of the matrices X^- and Xt and vectors {W\ and \V) appearing 
in (9.6-9.7) is proved by constructing them within an explicit representation. This 
representation has basis vectors that correspond to configurations of the rightmost 
N — k + 1 sites of the A'"-site lattice. These we denote as jT^Tfe+i • • -tjv). The vector 
\V) is then ascribed the role of a vacuum state, and the matrices Xr that of creation 
operators in such a way that 

Xr\V) = |t) and Xr\Tk+lTk+2 • • • Tjv) = \TTk+lTk+2 • • • Tjv) . (9.9) 

One then defines the vector {W\ via the scalar products 

{W\tiT2 ... tjv) = fin, T2, ■ • • , Tjv) (9.10) 
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so that (9.2) gives the desired statistical weights. 
The auxihary matrices Xt are defined as 

Xr ^ (^J^^kj+i + hn"^ X, (9-11) 

where the operators hjj+i and Iir are extended to the full space of all sub-configurations 
IrfcTfc + 1 • • • Tat) as follows. First, the bulk operator is defined in terms of the microscopic 
transition rates W{C ^ C) as 

(TfeTfc+i ■ ■ ■T-T-+^ ■ ■ ■ TN\h^+l\TkTk+l ■ ■ ■ TiTi+i ■ ■ ■ Tn) = W ijiTi+i ^'i^'i+l) (9-12) 

when i > k and (t/,t/^i) ^ (Ti,Ti+i), and 

{TkTk+l ■ ■ ■ nn+i ■ ■ ■ TN\hi^i+l\TkTk+l ■ ■ ■ TiTi+l ■ ■ ■ Tn) ^ - ^ W {TiTi+i ^ T-tI_^_i) (9.13) 

under the same condition that i > k. All other elements of these operators are zero. 
Similarly, at the right boundary we have 

(-Tfe • --T'l^lhRlrk • • - Tat) = W{tn ^ r'j^) (9.14) 

when Tn ^ t'n, and 

{Tk---rN\hR\Tk---TN) = - ^W{tn -^t'n) . (9.15) 

With these definitions established, the relations (9.3)-(9.5) follow after some 
straightforward manipulations [92]. 

Let us take stock of this construction. We have seen that not only does the existence 
of a set of vectors and matrices that satisfy (9.3)-(9.5) imply a stationary solution of 
matrix product form of the master equation (9.1), but also that such a set of vectors and 
matrices can always be constructed once one knows the weights (9.10). As we previously 
saw in Section 5 for the PASEP, one can find a number of different representations of 
the matrices D and E (which correspond to Xi and Xq in the more general setting) even 
when the auxiliary matrices D and E are the same. In the representation detailed above 
the auxiharies are not scalars but rather comphcated objects. Hence we see that (i) there 
may be many choices of both the matrices Xj- and their auxiliaries that correspond to 
the stationary solution of a single master equation; (ii) matrix reduction relations, like 
those found for the ASEP (2.39)-(2.41), constitute only sufficient conditions on X^., 
(W\ and \ V), since valid representations where these relations do not hold can be found 
(that used in this section provides an example); and (iii) the construction of the auxiliary 
matrices through the generators of the stochastic process in (9.11) does not necessarily 
imply the existence of any convenient reduction relations for the X matrices that allow, 
for example, efficient computation of statistical properties. 
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9.1.1. Existence of a matrix product state for models on the ring As we saw in Section 4 
the matrix product state on the ring involves using a trace operation. The algebraic 
relations to be satisfied are given in the general case by (9.6). However, as we saw in 
Section 7, sometimes it happens that, although the algebraic relations are consistent, 
the rotational invariance of the trace operation leads to global constraints on particle 
numbers. For example, in the ABC model discussed in 7.1.4 we found that the matrix 
product state was consistent only for models in which the numbers of each particle 
species was the same. It is essentially this property of models with periodic boundary 
conditions which has precluded the development of an existence proof for a matrix 
product state parallel to that given above. These difficulties are discussed in more 
detail in [168]. 



9.2. N on- conserving models with finite- dimensional representations 



We now discuss a specific reaction-diffusion model whose steady state can be represented 
in matrix product form, but where the auxiliaries are matrices. This model contains 
one species of particle and the stochastic processes are diffusion, coagulation and 
dccoagulation. These last two updates involve the annihilation or creation of a particle 
adjacent to another particle. If the bulk rates are chosen in the following way [169] 



11^01 

Aq 



01^10 11^10 

then for suitable boundary conditions the steady state may be written in matrix product 
form. 

Originally in [169] a closed segment was considered. Then, one requires D\V) — 
E\V) — (W\D — (W\E — the conditions to be satisfied come from the bulk algebra 

EE -EE =0 

ED-ED = q-^DE + q'^DD - (A + l)qED 
DE -DE = qED + qDD - (A + l)q-^ED 
DD- DD = AqED + Aq'^DD - {q + q'^)ED 

A four-dimensional representation of the matrices was found 
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q-^ 
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A/(A + 1) A/(A + 1) 
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with boundary vectors 

{W\ = (1 - q\ 1, 0, a) \V) = {b, 0, g^ - 1)^ (9.16) 

where a ^ b so that (W^IV) ^ 0. 

However there is a subtlety with this model on a closed segment which is that the 
empty lattice is dynamically inaccessible and itself comprises an inactive steady state. 
If we wish to exclude this configuration we should choose a and b so that 
vanishes. This can be accomphshed but in doing so renders a and b A^-dependent. 

A first order phase transition occurs at g^ = 1 + A. For g^ < 1 + A the system is 
in a low-density phase whereas for > 1 + A the system is in a high-density phase. 
Within the matrix product form of the steady state the transition can be understood 
as the crossing of eigenvalues of C. This may occur here not because the matrices are 
of infinite dimension but because they have negative elements (cf. the transfer matrices 
that arise in equilibrium statistical mechanics, as described in Section 2.3.3). 

Various other non-conserving models with finite-dimensional matrix product states 
have been identified. For example, the bulk dynamics just described has also been 
studied in the case of an open left boundary, where particles enter with rate a and 
leave with rate 7 and a closed right boundary by Jafarpour [170]. Then, one requires 
D\V) = E\V) = and 

-{W\D ^ a{W\E - -f{W\D (9.17) 

-{W\E ^ -f{W\E - a{W\D . (9.18) 

Then if 

a=(^-q + p^A (9.19) 

one can find two-dimensional representations of the algebra. Further generalisations 
of this class of model have been shown to have finite-dimensional matrix product 
states [171]. A general systematic approach to determine a necessary condition for the 
existence of a finite-dimensional matrix product state has been put forward by Hieida 
and Sasamoto [172] in which more examples are given. 

9.3. Several classes of particles 

Second-class particles were discussed in Section 2.4 and the matrix product solution to 
the second-class particle problem was detailed in Section 4. The matrix algebras for two- 
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species generalisations of the second-class particle arose in the BI and BII solution classes 
of the classification scheme of Section 7. We now consider a natural generalisation of the 
second-class particle to many classes of particle, labelled t — 1,2, . . . ,n with dynamics 

tO^Ot Vt and tt'^t't if r<T'. (9.20) 

In this case it turns out one requires more complicated operators and auxiliaries than 
matrices. 

For n — 3, the case of three particle classes, the algebraic relations to be satisfied, 
(9.6), become 

XiXi-XiXi^O, i = 0,l,...3 (9.21) 

and 



Xt-Xq 

X1X2 
X1X3 



XoX-r 
X2X1 
X3X1 
X3X2 



XqXt- 

X2XI 
X3X1 
^3X2 



-XrXo + XrXo , r > 
-XIX2 + XIX2 
-X1X3 + X1X3 
-X2X3 — X2X3 . 



(9.22) 



The two-class problem, where the last lines of (9.22) are absent, can be solved (as we 
have seen in Section 4.2) with scalar choices Xi = —1, Xq = 1 and X2 = 0. It turns out 
that three-class problem is not a simple generalisation of the two-class case, but instead 
a much more complication solution is needed [173]. Matrices which satisfy (9.21,9.22) 
have been given as 
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(9.23) 



The objects D and E that appear as elements of these infinite-dimensional matrices are 
themselves infinite dimensional matrices D and E which satisfy the familiar relation 
DE = D + E. In other words the Xi are rank four tensors, as indeed are the auxiliaries 
in this representation. These were also given by [173], and take the form 



/ D/2 + 1 E/2 - 1 
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(9.24) 



The tensor, or matrix within matrix, form of (9.23) is closely related to that of the 
operators used to determine the steady state current fluctuations in an open system 
[164]. 

In [173], some algebraic relations involving only the matrices Xj. (i.e., not the 
auxiliaries) were quoted. These included relations which had triples and quartets of 
matrices reducing to pairs, and others which transformed triples and pairs to other 
triples and pairs. As noted in [173] , the relations given comprise only a subset of those 
that would actually be required to perform a complete reduction of any string of matrices 
to an irreducible form. Hence, a concise statement of the steady state of this system, 
like that encoded by the reduction relations (2.39)-(2.41) for the ASEP, is not easy. 
As such, no calculations beyond simple cases [173] have yet been attempted for the 
stationary properties of this model. 

However recently, generalising a construction by Angel [64] for the two-class 
problem, Ferrari and Martin [174] have shown how to generate the steady state of 
the multi-class problem. This construction can then be used to determine operators 
and auxiliaries to solve the steady state for more than three classes [175]. 



10. Discrete-time updating schemes 

So far we have reviewed interacting particle systems involving continuous time, or 
equivalently, random sequential dynamics as discussed in Section 2. However, this is 
not necessarily the most natural choice of dynamics with which to model some physical 
systems of interest. For example, in traffic flow and pedestrian modelling [33, 50] it is 
desirable that the microscopic constituents are able to move simultaneously — this often 
originates from the existence of a smallest relevant timescale, e.g. reaction times in 
traffic. Therefore, in these systems, parallel dynamics in which all particles are updated 
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Figure 30. Typical configuration in the steady state of the partiaUy asymmetric 
exclusion process on a closed segment of TV sites with q < 1. 



in one discrete timestep, may be used. In this section we consider the ASEP under three 
types of discrete-time dynamics: sublattice, ordered sequential and fully parallel. 

In all these discrete-time cases the steady state weights are given by the eigenvector 
with eigenvalue one of the transfer matrix of the dynamics, which we write in a schematic 
notation is 

f/(ri,...,r;v) = /(ri,...,r^). (10.1) 

The transfer matrix T applied to the weight /(ti, . . . , r/v) generates a sum of terms, 
each being the weight of a configuration multiplied by the transition probability in one 
discrete timestep from that configuration to {ri, . . . , tn}. 



10.1. Sublattice parallel dynamics 

In sublattice updating the timestep is split into two halves: in the first half site 1, 
the even bonds (2, 3), (4, 5) . . . (2L — 2, 2L — 1) and site 2L are updated simultaneously; 
in the second half time step the odd bonds (1, 2),(3, 4) . . . (2L — 1,2L) are updated 
simultaneously. Note that we require the total number of sites N = 2L to be even. 

In Figure 30 we show how this type of updating is applied to partially asymmetric 
exclusion dynamics with open boundary conditions. In an update of site 1, if site 1 is 
empty a particle enters with probability a. In an update of a bond i,i + 1, the two 
possible transitions are: if site i is occupied and site i -|- 1 is empty the particle moves 
forward with probability p, or else if site i + 1 is occupied and site i is empty the particle 
moves backward with probability q. In an update of site A^, if site is occupied the 
particle leaves the system with probability /?. We note that this two-step process avoids 
the possibility of a confiict occurring (e.g., two particles attempting to hop into the same 
site simultaneously). 

The ASEP with sublattice dynamics was first studied in the case of deterministic 
bulk dynamics (p = 1, g = 0) by Schiitz [176]. A matrix product solution for this case 
was found by Hinrichsen [177] and a matrix product solution for the general case of 
stochastic bulk dynamics was found by Rajewsky et al [76,163]. We outline this general 
case here. 

One proceeds by constructing the transfer matrix for the whole timestep as a 
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product of operators for each half-timestep 



(10.2) 



where 



Ti 
T2 



LT23T45 
TuTu ■ ■ 



■ ■T2L-22L-1R (10.3) 

■f2L-12L (10.4) 

Tjj+1, L, R when acting on /(xi . . . r^r) generate the weights of configurations from which 
Ti . . .tn could have been reached by a transition associated with that bond or boundary 
site, multiphed by the transition probabihty. 
We may write 
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(10.6) 



where the basis for f is 00,01,10,11 and for L and R is 0,1. So, for example, 

f;,+i/(---T,_iiOTi+i---)= (10.7) 

pf{- ■ ■ Ti_i 1 Ti+i ■■■) + {!- p)fi- ■ ■ Ti_i 1 Ti+i ■■■). 

The matrix product state for sublattice updating takes the form 

/iv(ri, T2,..., Tn) = {W\Xr,Xr, ■ ■■Xr,,_,XrJV) (10.8) 

in which, for i odd, X^^ is a matrix D if site i is occupied by a particle (rj = 1) or a 
matrix E otherwise; and for, i even, Xri is a matrix D if site i is occupied by a particle 
or a matrix E otherwise. Note that different matrices, hatted and unhatted, are used 
for the odd and even sublattices. 

The cancellation mechanism is as follows 



-ii+l 



X-t-Xt- 



Xt-.Xt-. I -, 

RX^ 



(10.9) 

{w\LXr. = mxr, RXr^AV) = KAV)- (10.10) 

Note that the action of the first half time step transfer matrix Ti is to put hatted 
matrices at the even sites and unhatted matrices on the odd sites, then the action of T2 
is to restore them to their original sites. Thus, the matrix product state (10.8) is indeed 
an eigenvector with eigenvalue one of the full transfer matrix (10.1). Using the form 
(10.5,10.6) of the operators, the mechanism (10.9, 10.10) imphes the following algebraic 
relations 

[E,E] = [D,b] = 

ED = (1 - q)EID + pDE (10.11) 
IDE = {l-p)DE + qED 
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and the boundary conditions 

{W\E{l-a) = {W\E {l-m\V) = D\V) 

{W\{aE + b) = {W\D {E + (3D)\V) = E\V) 

Remarkably, the same matrices and vectors D, as in the random sequential case can 
be used to solve these relations. To see this we make the ansatz [163] 

E = E + \\ ■ D = D-X1, (10.13) 
where A is a scalar to be fixed, and we find that (10.11, 10.12) reduce to 

pDE-qED^X[{l-p)D + {l-q)E] (10.14) 
D\V) = ^\V) ; {W\E ^ {W\^^^^. (10.15) 
Finally we may set 

and define 

E'^X{l-q)E ; D' = X{l-p)D (10.17) 
which then satisfy 

pD'E' - qE'D' ^D' + E' (10.18) 

Thus, we have rewritten (10.14,10.15) in terms of matrices obeying the usual algebra for 
the PASEP (5.1)-(5.3) with redefined a and /3. From this rewriting, the phase diagram 
can, in principle, be deduced using the general results of [115]. 

10.2. Ordered sequential dynamics 

Another discrete-time updating scheme is to update each site in a fixed sequence in each 
time step. Two particularly obvious choices of sequence are as follows. 

Forward- ordered update Here, the timestep begins by updating site 1 wherein if site 
1 is empty a particle enters with probability a. Next, the bond 1,2 is updated such 
that if site 1 is occupied and site 2 is empty the particle moves forward with probability 
p, or else if site 2 is occupied and site 1 is empty the particle moves backward with 
probability q. Then the bonds i, i + 1 ior i — 2 . . . N — 1 are similarly updated in order. 
The time step concludes with site N being updated wherein if site is occupied the 
particle leaves the system with probability f3. Note that in the forward sequence it is 
possible for a particle to move several steps forward in one timestep. 
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Backward-ordered update In this case, the time step begins by updating site N, then 
bonds iV — 1, to 1, 2 in backward sequence and finally site 1, where all the individual 
updates follow the usual rules, as described above. Note that due to particle-hole 
symmetry the dynamics of the vacancies in the backward order is the same as the 
updating of particles in the forward order. Therefore there is a symmetry between the 
steady states for the forward and backward orders. 

Let us focus on the backward-ordered dynamics. We construct the transfer matrix 
for the whole timestep as a product of operators for each update of the sequence 

f = Lfi2f23 . . . fW-lAT^ (10.20) 

where L, Tjj+i, R are as in (10.5,10.6). 

The matrix product solution is of the usual form 

/iv(ri, r,,..., Tn) = {W\X,,X,, ■■■X^^\V) (10.21) 

The cancellation mechanism is precisely the same as for the sublattice parallel case 
(10.9, 10.10). In the ordered case the hat matrices do not appear in the steady state 
weights at the end of the time step, rather they are auxiliary matrices which appear 
during the update procedure and move from right to left through the lattice. 

We conclude that backward ordered updating has the exact same phase diagram 
as sublattice updating, although exact expressions for correlation functions do depend 
on the details of the updating [76, 178]. Finally the phase diagram for forward ordered 
updating can be obtained from the particle-hole symmetry mentioned above. Current 
and density profiles have been calculated in [179,180]. 

10.3. Fully parallel dynamics 

In parallel dynamics (sometimes referred to as fully parallel to distinguish from the 
sublattice updating) all bonds and boundary sites arc updated simultaneously. This 
dynamics is considered the most natural for modelling traffic fiow [33]. In an update 
at the bond + 1), if site i is occupied and site i + 1 is empty the particle moves 
forward with probability p. Under this type of updating scheme one cannot include 
reverse hopping without introducing the possibility of conflicts occuring, or of a particle 
hopping to two places at once. Also note that it is the occupancies at the beginning of 
the timestep which determine the dynamical events. 

In this case the matrix product solution is of the usual form (10.21) but the algebraic 
relations have a more complicated structure which was first elucidated in [178]. There 
it was found that recursion relations between systems of different sizes were higher than 
first order: for example the relation 

/jv(. . . 0100 . . .) - (1 -p)/jv-i(. . . 010 . . .) + /jv-i(. . . 000 . . .) +pfN-2{. . . 00 . . .) ,(10.22) 

which relates the weights of size to those of size — 1 and size N — 2, was discovered 
to hold. This in turn implies that algebraic relations between the operators are quartic 
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rather than quadratic. For example the rules in the bulk are 

EDEE ^{l-p)EDE + EEE+pEE (10.23) 

EDED = EDD + EED + pED (10.24) 

DDEE = (1 - p)DDE + (1 - p)DEE + p(l - p)DE (10.25) 

DDED = DDD + {1 - p)DED +pDD (10.26) 

and there are other rules which we do not quote here for matrices near to the boundary 
(see [178]). These rules were proved using a domain approach related to, but more 
complicated than, that presented in in Section 3.1 [178]. Subsequently, an algebraic 
proof similar in spirit to the cancellation mechanism for the ordered sequential updating 
case was found [181] but it is still too involved to present here. 

The quartic algebra, (10.23)-(10.26), as well as the other conditions mentioned 
above, can be reduced to quadratic rules by making a convenient choice for the operators 
involved. The trick is to write 



r>=r' : ^ E= „■ ? . (10.27) 





where Di, D2, Ei, and E2 are matrices of, in general, infinite dimension; that is, D and 
E are written as rank four tensors with two indices of (possibly) infinite dimension and 
the other two indices of dimension two. Correspondingly, we write {W\ and \V) in the 
form 

{W\^{{W,\, {W,\), \W)^(^ jj^^l j , (10.28) 

where {Wi\, (W2I, |Vi), and IV2) are vectors of the same dimension as Di and Ei. Di, 
El, {Wi\, and |Vi) satisfy the quadratic relations 

DiEi = (1 - p) [Di + Ei+p] , (10.29) 

^il^i) = ^^^^^V i) , {Wi\E, = {W,\ ^^^~''\ (10.30) 
p a 

and D2, E2, (VF2I, and |V^2) satisfy 

E2D2^p[Di + Ei+p], (10.31) 

E2\V2)=p\Vi) , {W2\D2 = {Wi\p. (10.32) 

and {Wi\ satisfying (10.29,10.30) are presented in [178]. In addition, one can choose 
D2 oc El, E2 oc Di, {W2\ oc (H^il, 11/2) oc \Vi) so that (10.31), (10.32) reduce to (10.29), 
(10.30). Along the curve 

(1-«)(1 -/?) = 1 -p, (10.33) 

scalar representations of Di, Ei, \Vi) and {Wi\ can be found and D, E are 2x2 operators. 

The continuous time limit is recovered by setting p = dt, replacing a — > adt, 
j3 j3dt and letting the time step dt — >■ 0. Then (10.29,10.30) reduce to the usual 
ASEP quadratic algebra (2.39)-(2.41). 
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Figure 31. Phase diagram for the ASEP with parallel update. Note that a < 1, 
/3 < 1. MC is the maximum current phase, LD and HD are the low and high-density 
phases, respectively. The straight dashed lines are the boundaries between the low 
density and maximal current phases and the high density and maximal current phase 
at a = 1 — y/1 — p and /? = 1 — y^l — p respectively. The curved dashed line is the line 
given by (10.33) and intersects the line a = /3ata = /3=l — y/l — p = q 



The phase diagram is presented in Figure 10.3. It appears similar to the famihar 
continuous time phase diagram except that the transition hues from low density to 
maximal current and high density to maximal current are at a = 1 — y/1 — p and 
j3 = 1 — a/I — p respectively. In the continuous time limit described above one recovers 
a = 1/2 and f3 = 1/2. Also, in the limit of deterministic bulk dynamics, p — 1 [182] the 
maximal current phase disappears from the phase diagram as is also the case in other 
updating schemes in this limit. 

11. Summary and outstanding challenges 

In this work we have reviewed the physical properties of systems that are driven out 
of equilibrium, with particular reference to those one-dimensional models that can be 
solved exactly using a matrix product method. The most prominent physical phenomena 
exhibited by these systems are phase transitions which in open systems are induced by 
changes in the interactions with the boundaries, and in periodic systems by the addition 
of particle species that have different dynamics. Additionally we have seen in a number 
of cases the formation of shocks in these one- dimensional systems. 

As we have also discussed, the matrix-product approach can be used to determine 
the steady-state statistics of models with a variety of microscopic dynamics. The 
majority of models that have matrix product states involve diffusing particles with 
hard-core interactions and whose numbers are conserved, except possibly at boundary 
sites. Under such conditions, certain generalisations to multi-species models have been 
found. Furthermore, there are a few cases in which non- conserving particle reactions 
can occur that also admit a convenient solution in terms of matrix products. 

However, it has to be said that the dream scenario, of being able to systematically 



108 



construct any nonequilibrium steady state in matrix product form starting from the 
microscopic dynamics of the system, appears a remote goal. Ahhough the existence 
proofs discussed in Section 9 teU us that a matrix product formulation should nearly 
always be possible, we are in practice still feeling our way with particular examples. 

There remain a number of simple physical systems that exhibit nontrivial out-of- 
equilibrium behaviour, but that have nevertheless been unyielding to exact solution by 
matrix products, or any other means. It is perhaps appropriate here to highlight some 
of these cases as challenges for future research. 

11.1. Particle- and site-wise disorder in the ASEP 

As was discussed in Section 8.1, it is possible to solve for the steady state of the 
ASEP on the ring when particles have different hop rates, as would occur when each 
particle is randomly assigned a hop rate from some distribution. As far as we aware, 
the corresponding dynamics with open boundaries has not been solved apart from the 
Karimipour's model of overtaking which we reviewed in Section 8.2. Furthermore, if the 
disordered hop rates are associated with sites, rather than particles, even the model on 
the ring has so far eluded a complete exact solution. Indeed, even if only a single site has 
a different hop rate to the rest, the steady state does not appear to have a manageable 
form [183, 184] . On the other hand, some nontrivial symmetries in the disordered case 
have been identified [185, 186]. 

One phenomenon that can emerge when disorder is present is phase separation. 
One sees this most clearly for the periodic system with particle-wise disorder, where as 
we noted in Section 8.1.1, a condensation of particles behind the slowest may occur [140]. 
With site- wise disorder [185], a flattening of the macroscopic current-density relation 
J(p) is seen (recall that it is parabolic J(p) = p(l — p)in the ASEP without disorder), 
and the densities at which transitions into a maximal current phase are correspondingly 
lowered [187]; it has been suggested that a linear portion in the J(p) curve may also 
provide a mechanism for phase separation [188] . Another interesting observation is that 
under site-wise disorder the location of the first-order transition for the ASEP becomes 
sample-dependent (i.e., the free-energy-like quantity for the nonequilibrium system is 
not self-averaging) [189]. 

11.2. Driven n-mers with open boundaries 

In all the models described in this work, particles occupy a single site of a lattice. If one 
has extended objects on a ring, the dynamics are unaffected. With open boundaries, 
however, the situation changes and even if there is a single particle species and all 
particles are the same length, a matrix product solution for the steady state has proved 
elusive. Although it is quite easy to arrange for cancellation of terms in the master 
equation corresponding to particles joining and leaving domains (see Section 3.1), the 
necessary cancellation with boundary terms is much harder to arrange: it is not even 
clear what the most convenient choice of boundary conditions should be (e.g., extended 
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particles sliding on and off site-by-site, or appearing and disappearing in their entirety, 
or a mixture of the two). A phcnomcnological diffusion equation for the macroscopic 
density profile in the open system with extended objects appearing and disappearing 
at both boundaries has been presented [190]. Furthermore, it was conjectured in that 
work to be the equation that would result in the continuum limit (along the lines of 
that taken in Section 2.2), if an exact matrix-product solution of this model were to 
be found. Whether this should turn out to be the case or not, an apphcation of the 
extremal current principle (see 2.2.4) gives predictions for the phase diagram that are in 
good agreement with simulation data: the same three phases that emerge in the ASEP 
are seen, but with shifted transition points. Meanwhile an approach to this problem 
based on a local-equilibrium approximation [191] has shown excellent agreement with 
Monte Carlo results and have also been conjectured to be exact. The case of the open 
system with extended objects and a single spatial inhomogeneity in the hop rates has 
also been considered [192]. 

11.3. Bridge model 

The bridge model [145, 146] comprises oppositely-charged particles in the sense of 
Section 7. Positive charges enter at the left at rate a, hop with unit rate to the right 
and exit at the right boundary at rate (3. The dynamics of negative charges is exactly 
the same, but in the opposite direction. When two opposing particles meet they may 
exchange with some rate q. Thus this model has both a parity and charge-conjugation 
symmetry. In certain parameter ranges, this symmetry is reflected in the steady state 
that results: both positive and negative charges flow freely in their preferred directions. 
In particular, the limit a ^ oo at j3 = 1, for which an exact solution in terms of 
matrix products exists, is included within this range. On the other hand, when f3 is 
small the model's symmetries arc found to be broken. This has been proved in the 
(3^0 limit [193] but there is no known exact solution for general f3. In a finite system, 
the current flips between the positive and negative direction, but as the system size 
increases the time between flips increases exponentially, and so in the thermodynamic 
limit only the state with positive or negative current state is seen. It is unclear whether 
the structure of the matrix product solution allows for the description of such symmetry- 
broken states. This may also be related to whether the boundary conditions correspond 
to particle reservoirs with well-deflned densities [194] 

11.4. ABC model 

The ABC model [141, 142] is another deceptively simple model that has yet to be 
completely solved. This model has a ring that is fully occupied by particles from three 
species, which are labelled A, B and C and exhibit a symmetry under cyclic permutation 
of the labels. Specifically, the rates at which neighbouring particles exchange arc 

AB^BA BC^CB CA^AC. (11.1) 
111 ^ ' 
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That is, if g < 1, A particles prefer to be to the left of Bs. Bs to the left of Cs and Cs to 
the left of As. This model was discussed in Section 7.1.4 and was shown to be exactly 
solvable when the numbers of A, B and C particles are all equal as a consequence of 
detailed balance being satisfied. Here a phase separation is seen when q ^ 1. If g = 1, 
on the other hand, all configurations become equivalent and the stationary distribution 
is uniform. When there is even the slightest imbalance and particle numbers and hop 
rates, a small current starts to flow and nothing is known about the properties of the 
stationary distribution. 

The model is of particular importance because in the special case of equal numbers 
of each species, it is one of few models along with the SSEP and KMP model [195] where 
a free energy functional for the density profile is known. 

11.5. KLS model 

We end by mentioning the original paradigm of a system driven out of eqiiilibrium, 
the model due to Katz, Lebowitz and Spohn [30,196,197]. This model has Ising-like 
interactions between pairs of neighbouring spins and evolves by neighbouring spins 
exchanging places (Kawasaki dynamics) but with a symmetry-breaking external field 
that favours hops of up-spins along a particular axis. When the system is open, or has 
periodic boundaries, a nonequilibrium steady state ensues. In the limit where spins 
become non-interacting and there is only one spatial direction, the partially asymmetric 
exclusion process (PASEP) described in Section 5 is recovered. When the spins interact, 
an Ising-like steady state that involves products of 2 x 2 matrices (as described in Section 
2.3.3) can be found, albeit only if certain conditions on the hop rates are satisfied. (These 
conditions do admit nonequilibrium steady states in which a current fiows, however). 
In two dimensions, one has, of course, a finite critical temperature in the absence of a 
driving force, and there has been particular interest in the nature of the low-temperature 
ordered phase in the presence of an external drive [30] . 

11.6. Fixed random sequence and shuffled dynamics 

We discussed in Section 10 different updating schemes for which matrix product steady 
states have been determined. Here we mention some other updating schemes which have 
not been solved so far. First an ordered discrete time sequential updating scheme could 
have the order as a randomly chosen sequence. The forward and backward schemes 
of Section refdiscretetime correspond to special sequences. It would be of interest to 
construct the a matrix product steady state for an arbitrary fixed sequence. 

Recently a shuffled updating scheme has been considered [198, 199], in which a 
new random sequence is chosen at each timestcp. It has been argued that this scheme 
is relevant to the modelling of pedestrian dynamics. The shuffled updating scheme 
guarantees that each site or bond is updated exactly once in each timestcp but is more 
stochastic than an ordered sequential scheme. Again it would be of interest to see if a 
matrix product could be used to describe the steady state. 



Ill 



Acknowledgements 



We thank all the many colleagues with whom it has been a pleasure to discuss 
nonequilibrium matrix product states over the years. In particular we acknowledge 
our collaborators and other authors whose work we have summarised. For their useful 
discussions and insightful comments during the preparation of this manuscript we 
would wish to thank Bernard Derrida, Rosemary Harris, Des Johnston, Joachim Krug, 
Kirone Mallick, Andreas Schadschneider, Gunther Schiitz and Robin Stinchcombe. RAB 
further acknowledges the Royal Society of Edinburgh for the award of a Research 
Fellowship. 



A. Method of characteristics 

The continuity equation (2.3) is a particularly simple case of a first order quasi-linear 
differential equation for the density 

The term on the right hand side would represent a source or sink of particles and is 
relevant for example to the case where creation and annihilation processes exist. Such 
an equation can generally be solved by the method of characteristics [200] which involves 
identifying the characteristic curves along which information from the boundary or initial 
condition propagates through the space-time domain. The characteristic curves satisfy 
dt dx dp / . X 

- = -r = — A.2 

a c 

of which the two independent conditions may be written 
da; ^ 6 dp ^ c 

dt a dt b ^ ' ^ 

These equations generally have two famihes of solution 

=Ci i;{x,t,p{x,t))^C2 (A.4) 

Then the solution to (A.l) is of the form 

F(0,V^) = O (A.5) 

where the function F is fixed by the initial data. 

For the case (2.3), i.e. c = 0, the characteristics in the x-t plane are straight lines 
with slope Vg{p) and p is constant along these fines, as can be seen from (A. 3). More 
generally the characteristics are curves in the x-t plane and the density varies along the 
characteristic. This can be seen by integrating (A. 3) 

G{p)-G{po) = t where G{p) = T dp^^ (A.6) 

which gives p implicitly as a function of t, then the trajectory of the characteristic is 
given by 

x^Xo+ [ bipit'))dt' . (A.7) 
Jo 
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If we consider a patch of the initial density profile of density po at Xo, its density evolves 
according to (A. 6) and (A. 7) gives the position of the patch. These more complicated 
characteristics can produce a richer for boundary induced phase transitions [43, 201] 
including stationary shocks. 



B. Generating functions and asymptotics 

In this appendix, we collect together general formulae for obtaining both exact and 
asymptotic expressions for the coefficients appearing in a generating function if the 
latter is known. We refer the reader to [95] for proofs and more detailed discussion. 
For brevity, we introduce the notation {x"'}f{x) to mean "the coefficient of x"' in the 
(formal) power series f{x)". 

Lagrange Inversion Formula If a generating function f{x) satisfies the functional 
relation 

f{x)=x^f) (B.l) 

where $(/) satisfies $(0) = 1, one has an expression for a; as a function of / which can 
be inverted to find / as a function of x. This procedure then allows one to determine 
with relative ease the coefficient of in the expansion of any arbitrary function F{f). 
The formula for doing so reads 

{x-}F{f) = -ir-'} [F'imfr] . (b.2) 

Th 

Asymptotics If a generating function f{x) has a singularity at a point Xq in the complex 
plane, it contributes a term Axq'^u'^ to the coefficient {x'^}f{x), where A and fi are 
constants to be determined. For sufficiently large n, this coefficient is dominated by the 
singularity closest to the origin. To be more precise, we say that two sequences a„ and 
bn are asymptotically equivalent, denoted Un bn, if 

lim ^ = 1 . (B.3) 

n^oo bn 

Then, we have that {a;"}/(a;) ~ Ax^'^n''^. 

Let us decompose the generating function f{x) into its regular and singular parts 
around Xq, denoting the former by fr{x). That is, 

f{x) = Mx) (l - ^) ' ■ (B.4) 

The cases that interest us are poles u — 1,2, . . . and algebraic singularities (noninteger 
v). For the former case one finds 

{x-}fix) ~ ^ ^ " ^) /.(^o)% " . (B.5) 

For sufficiently large n, the binomial coefficient behaves as 
V + n - 1\ n"-^ 
n J (i/ — 1)! 
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Replacing the factorial with a Gamma function gives the extension of (B.5) to algebraic 
singularities 

The pref actor fr{xo) in this expression can be obtained by taking the limit 

fr{xo)^ lim (l--) fix). (B.8) 

C. Equivalence of the integral and sum representations of the ASEP 
nor malisat ion 

In order to show that the integral representation of the normalisation (3.49) is equivalent 
to the summation (3.38), the simplest approach is to show that the generating functions 
for computed from the two expressions are equal. As we have seen in Section 3.2.3 
the generating function computed from the finite sum (3.38) yields (3.56) which can be 
developed further to give 

oo 

Z{z) = J2 ^""Z^ 

^ aP [2a(3 - 2z - (a + (3 - + VT^)] 

2[z-a{l-a)][z- ^)] ' ^ ' ^ 

Wc now turn to the integral representation (3.49) which may be recast as a contour 
integral by change of variable u = c*^ 

(1 - ab) f du (2 - _ u-^){2 + U + m"^)^ 



(C.3) 



Airi Jj^ u [1 — au){l — au — hu){l — hu ^) 
where K is the (positively oriented) unit circle in the complex u plane. 

We compute the generating function (C.l) using (C.3). Summing a geometric series 
(which converges for \z\ < 1/4) we find 

Z(z\ = - / ^ {2-v?- u-') 

^' Am Jk u {l-au){l-au-^){l-bu){l-bu-^){l- z{2 + u + u-^)) ^ ' ^ 

We now use the residue theorem to evaluate this integral. We take a and 6 < 1 (although 
we get the same final result for other ranges of a and b). In this case the singularities 
within the unit circle are simple poles at u — a, u — b and u — where we define 

«± = 1-^(1±VT^) (C.5) 
The residues from the three poles yield 

Z(z) = - - (^-^^) (c 6) 

2{a - b) [z - a{l - a)] 2{a - b) [z - - (3)] ^ ^ 

il-ab)i2-ul-u-J) 



2z{\ + a2 - a(M_ + m:^))(1 + 6^ - b[u- + mI^))(m- - «+) 
Some simple algebra then shows that this expression is equivalent to (C.2). 
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D. Matrix representations 



Here, we collect together some representations of the various matrix product algebras 
that have been discussed in this work. 



D.l. Representations of PASEP algebra 



We begin by recapping representations for the PASEP algebra in the case of injection 
at the left boundary and extraction at the right boundary 

DE-qED = D + E (D.l) 

(5D\V) = \V) (D.2) 

a{W\E - {W\ . (D.3) 

These representations are generalisations of the three representations for the totally 
asymmetric case first given in [58]. 



Prom Section 5.2 we have 

f 1 + b 

1 



D 



1-q 









V 



/ l + a 



E = 



1-q 








1 + bq 





1 + aq 
Vcl 






l + bq^ 







1 + aq'^ 






1 + bq^ 





1 + aq^ 



\ 



(D.4) 



(D.5) 



where 



1-q 

a 



- 1 



- 1 



l-g"+^)(l-a6g'^) 



and the boundary vectors are 

{W\ = (0| and \V) = |0) . 
From this representation one sees that for certain parameter curves, namely 
1 - a6g" = , 



(D.6) 

(D.7) 
(D.8) 



that the representations become finite dimensional since Cn — and the upper left corner 
of the matrices become disconnected from the rest. Curves of this type were first noted 
in the more general case ^,5 ^ in [202] and finite dimensional representations were 
catalogued in [203]. 



115 



In the limit ? — > 1 (5.12) and (5.12) have well-defined limits 
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1 + l/B 









D = 








2 + 1/3 


V • -' 















3 + 1//? 






V ; 








■■ / 




/ l/a 











• \ 






1 + l/a 










E = 







2 + l/a 

















3 + l/a 














■ / 



1 1 

where g — — \- —. 

a p 

A second representation of (D.1-D.3) given in Section 5.3 is 



/l 






1 















1 





v 



E 



1-q 



1 










1 






1 







1 



(D.9) 



(D.IO) 



(D.ll) 



(D.12) 



with 



{W\n) = K 



and {n\V) — k- 



(D.13) 



V(?; (i)n ' ' ' \/{q;q)7 

where the parameters a and h arc given in (D.6), and /t is a constant usually chosen so 
that = 1. For these matrices, the g — > 1 limit is not defined. 

A third representation, which we have not yet encountered, is 

/I//? 1//3 1//5 1//5 •••\ 

^(13) 
D = d(22) ^(23) 





v 



/ 
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/ ••• \ 

10 
E = 10 
10 

V; 

{W\n) = {l/aY {n\V) = 6^,0 



where 



i-l 



m=0 



j -i + m\ ^(t_n 
m J f3 \i 



In the q ^ 1 hmit the elements of D reduce to 



J 



(D.14) 

(D.15) 
(D.16) 

(D.17) 



D.2. Representations of two-species and multispecies algebras 

We now give representations corresponding to the physicaUy relevant and nontrivial 
algebras classified in Section 7 and the multispecies generalisations of Section 8. 

Solution All and multispecies generalisation This is a multispecies model with matrices 
E and D{v), with an algebra 



D{v)E ^-D{v) + E 



V 



D{v)D{v') = -D{v') 



D{v)\V) 
{W\E 



V — V 
1 



V — V 



-D{v) V > v' 



1 

a 



\V) 



(D.18) 



where j3{y) 



D(v) 



One possible representation is 
/l/(3{v) l/{P{v)v) l/{f3{vy) l/{f3{vy) 









1 






1/v 
1 




1/v 
1 



\ 



) 



E 



\ 



/ 
10 
10 
10 

V : ■■• / 

= (1/a)" {n\V) = 6^,0 
This is the generalisation of the third PASEP representation (D.15) in the case q = 0. 
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(D.19) 



Solution BI and BIT With the algebra (7.18), 
PiXiXo - QiXqXi = xoXi - xiXo 
P2X2X0 — q2XQX2 — X0X2 
P2X1X2 — q2X2Xi — — X1X2 

one has a representation 



where 



Xi= - 



Xi 



(P2 - 



Xo 
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■■ \ 
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■■ / 
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{Q2/P2f 
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{P2 - 52) 



(1 - {qi/PiT^') 
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■■ ^ 
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VC2 
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(D.20) 



(D.21) 



(D.22) 



(D.23) 



(D.24) 



As noted in the main text a representation for solution Bll is to use a presentation of 
the PASEP algebra and write X2 = \V){W\, i.e. a projector. 



Solution CII and multispecies generalisation This is the case of the ASEP with 
disordered hopping rates. The matrix algebra (8.2) 



can be represented via 



D, 



( dT^ 




dT' 


4°'^ 






























4''^ 



(D.25) 
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E 



where 



/ n 
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n 
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" \ 
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■■ 1 






ql 







(D.26) 



1+3 



(D.27) 



Solution D and Deformed commutators As noted in 7.2.2 the Solution class D operators 
can be constructed from tensor products of operators obeying deformed commutator 
relations. The deformed commutator is defined as 



XrX2 - rXaXi = 



(D.2^ 



where r is a deformation parameter: r = 1 corresponds to the usual commutator 
[Xi, X2] = 0. Deformed commutators appear in the analysis of the PASEP, and in some 
of the two-species models of Section 7. One possible representation of these matrices is 



Xi 







V : 




1 





J 



( 1 








r 


















\ 



(D.29) 



We note that both Xi and any product of Xx and X^ including at least one Xx are 
traceless. 
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